PRIMERA CLASE: REGRESIÓN LINEAL

Maestría en Investigación Operativa y Estadística

UTP

INTRODUCCIÓN

DOCENTES

Héctor Hernán Montes García.

Ingeniero industrial de la Universidad Tecnológica de Pereira.

Estudios en Maestría en Ciencias con orientación en Matemáticas

Desempeño laboral

Científico de datos con más de 3 años de experiencia en la construcción de modelos de aprendizaje automático y soluciones de minería de datos para obtener mejores conocimientos comerciales. Experiencia utilizando SQL, bibliotecas de Python (matplotlib, seaborn, flask, scikit-learn, pandas, numpy), modelos matemáticos y estadísticos para ofrecer soluciones robustas que agregan valor al negocio.

Ha trabajado para clientes del sector industrial y financiero en México y Colombia, en áreas tales como:

  • Mantenimiento predictivo
  • Diseño de campañas comerciales basadas en datos
  • Reconocimiento de imágenes
  • Modelos de procesamiento de lenguaje natural (NER).

Julián Piedrahíta Monroy

Ingeniero industrial.

Magíster en desarrollo agroindustrial

Universidad Tecnológica de Pereira

Desempeño laboral

  • Consultor en soluciones análiticas para instituciones educativas.
  • Analista de datos en el observatorio social de la UTP.
  • Actualmente diseñador de tableros informativos (dashboards) para la Universidad Pedagógica Nacional de Bogotá y la misma Universidad Tecnológica de Pereira.
  • Docente catedrático en el área de informática y estadística general.
  • Uso principal del lenguaje R y sus librerias Tidyverse, RMarkdown, RShiny y otras.

Contenido del curso

  • T1. Regresión Lineal, suposiciones y requisitos.
  • T2. Estimación de los parámetros e interpretación de los valores.
  • T3. Pruebas de hipótesis relacionadas con los parámetros.
  • T4. Evaluación de modelos de acuerdo a las suposiciones de normalidad e independencia.
  • T5. Transformaciones de las variables, y sus implicaciones en las pruebas de hipótesis.
  • T6. Series de Tiempo. Estacionariedad y cómo obtener una serie estacionaria a partir de una que no lo es. Correlogramas.
  • T7. Modelos de Box y Cox. Estimación de parámetros.
  • T8. Criterios de evaluación de un modelo de series de tiempo. Transformaciones.

Bibliografía.

Regression Modeling and Data Analysis with Applications in R. Samprit Chatterjee, Jeffrey S. Simonoff
Montgomery, D. C., Peck, E. A., & Vining, G. G. (2002). Introducción al análisis de regresión lineal (3ª ed.). Limusa Wiley.
Time Series Analysis and Its Applications: With R Examples by Shumway and Stoffer.
Time Series Analysis: Forecasting and Control by Author(s): George E. P. Box, Gwilym M. Jenkins, Gregory C. Reinsel, Greta M. Ljung.
https://fhernanb.github.io/libro_regresion/rls.html#modelo-estad%C3%ADstico
https://bookdown.org/victor_morales/SeriesdeTiempo/

Motivación

Los métodos de regresión son técnicas estadísticas utilizadas para modelar y comprender la relación entre variables. La motivación detrás de estos métodos radica en la necesidad de analizar y cuantificar cómo una o más variables predictoras influyen en una variable respuesta. Los métodos de regresión son ampliamente utilizados en diversos campos, incluyendo la investigación científica, la economía, la medicina, la ingeniería y más. A continuación, se detallan algunas de las principales motivaciones detrás de estos métodos:

Img 1: Ejemplo de fenómeno lineal

Modelado de Relaciones:

La motivación fundamental de la regresión es modelar la relación entre variables. En muchos casos, existe la necesidad de entender cómo una variable dependiente cambia en función de una o más variables independientes. Por ejemplo, en la economía, podría ser necesario entender cómo las tasas de interés afectan el gasto del consumidor.

Predicción y Estimación: Los modelos de regresión también se utilizan para hacer predicciones y estimaciones. Dado un conjunto de datos históricos, los modelos de regresión pueden ayudar a predecir valores futuros de la variable dependiente en función de los valores de las variables independientes. Esto es especialmente útil en campos como la meteorología, la finanzas y el marketing.

Img 2: Prediciendo la tendencia histórica del S&P

Control y Optimización: En algunos casos, se utilizan modelos de regresión para optimizar y controlar procesos. Por ejemplo, en la manufactura, los modelos de regresión pueden usarse para identificar las condiciones ideales que conducen a la máxima eficiencia o calidad del producto.

Validación de Teorías: En la investigación científica y social, los modelos de regresión pueden usarse para validar o refutar teorías existentes. Al comparar los resultados del modelo con las expectativas teóricas, se puede evaluar la validez de las hipótesis.

Gestión de Riesgos: Los modelos de regresión también se utilizan para evaluar riesgos y tomar decisiones informadas. En finanzas, por ejemplo, se pueden utilizar para evaluar el riesgo de inversión y la exposición a diferentes factores del mercado.

Ejemplo

El modelo CAPM (Capital Asset Pricing Model) es un modelo de valoración de activos financieros desarrollado por William Sharpe que permite estimar su rentabilidad esperada en función del riesgo sistemático. Su desarrollo está basado en diversas formulaciones de Harry Markowitz sobre la diversificación y la teoría moderna de portafolios. Dos puntos claves son:

  1. El modelo CAPM es utilizado para calcular la rentabilidad que un inversionista debe exigir al realizar una inversión en un activo financiero en función del riesgo que está asumiendo.
  2. El modelo CAPM establece una relación lineal entre el rendimiento esperado de un activo y su riesgo sistemático, medido por su \(\beta_i\).

La fórmula del modelo CAPM es la siguiente:

\(E(r_i)=r_f+\beta_i*(E(r_m)-r_f)\)

Donde:

\(E(r_i)\) es la tasa de rentabilidad esperada de un activo concreto. \(rf\) es la rentabilidad del activo sin riesgo. \(\beta_i\) es la medida de la sensibilidad del activo respecto a su benchmark. \(E(r_m)\) es la tasa de rentabilidad esperada del mercado en que cotiza el activo.


Img 3: Modelo CAMP para valoración de activos financieros

Análisis de Causa y Efecto: Los modelos de regresión pueden ayudar a establecer relaciones causa-efecto entre variables. Esto es útil para comprender cómo los cambios en una variable influyen en otras variables y viceversa.

Img 4: Predicción de rendimientos de cultivos

En resumen, la motivación detrás de los métodos de regresión radica en la necesidad de entender, modelar, predecir y cuantificar las relaciones entre variables en una variedad de campos. Estos métodos permiten tomar decisiones informadas, realizar análisis profundos y obtener información valiosa a partir de los datos disponibles.

1. INTRODUCCIÓN NO FORMAL AL MODELO LINEAL SIMPLE

En este capítulo introduciremos el modelo lineal simple partiendo de datos reales del precio del aguacate y el precio del dólar. Más que una presentación formal o matemática del modelo, nuestro objetivo será construirlo paso a paso usando las ideas estadísticas detrás del método. Confiamos en que esto le al lector una idea más precisa de cómo éste modelo usa los datos dados para relacionar las variables de inteŕes, y que comprenda sus fortalezas pero también sus limitaciones.

1.1 Lectura de datos

A continuación, cargaremos dos tablas de datos que contienen información sobre el precio del dolar, también conocida como tasa representativa del mercado, y otra con la información del precio del aguacate Hass tomada de una fuente gubernamental y de la página del SIPSA.

Los enlaces para acceder a información relacionada son:

https://www.agronet.gov.co/estadistica/Paginas/home.aspx?cod=11 https://www.bde.es/webbe/es/estadisticas/temas/tipos-cambio.html

# Cargue y limpieza de datos.

dolar <- read.csv2("datasets/Tasa_de_Cambio_Representativa_del__Mercado_-Historico.csv") %>% 
  clean_names()
hass <- read.xlsx("datasets/Hass_Precios_Historicos.xlsx") %>% clean_names()

1.2 Creación de funciones.

Se programa una función de flextable que mejorará la impresión de tablas de resumen. Para usarla sólo bastará invocar la función ftable() al final de la sentencia (usando %>% tidyverse) o con el ftable() encapsular o encerrar la tabla que se quiere ajustar.

# Funciones

## Función para crear flextable
ftable <- function(x) {
  x %>% 
    flextable() %>% 
    theme_vanilla() %>%
    color(part = "footer", color = "#666666") %>%
    color( part = "header", color = "#FFFFFF") %>%
    bg( part = "header", bg = "#2c7fb8") %>%
    fontsize(size = 11) %>%
    font(fontname = 'Calibri') %>%
    # Ajustes de ancho y tipo de alineación de las columnas
    set_table_properties(layout = "autofit") %>% 
    # width(j=1, width = 3) %>%
    align(i = NULL, j = c(2:ncol(x)), align = "right", part = "all")
}

Como ejemplo del uso de la función:

Sin flextable

dolar %>% 
  head(5) 
##      valor unidad vigenciadesde vigenciahasta
## 1 4,761.64    COP      22/12/22      22/12/22
## 2 4,781.28    COP      20/12/22      20/12/22
## 3 4,802.48    COP      17/12/22      19/12/22
## 4 4,836.24    COP      13/12/22      13/12/22
## 5 4,815.99    COP      10/12/22      12/12/22

Con flextable

dolar %>% 
  head(5) %>% 
  ftable()

valor

unidad

vigenciadesde

vigenciahasta

4,761.64

COP

22/12/22

22/12/22

4,781.28

COP

20/12/22

20/12/22

4,802.48

COP

17/12/22

19/12/22

4,836.24

COP

13/12/22

13/12/22

4,815.99

COP

10/12/22

12/12/22

1.3 Depuración de los datos

Dolar

Como en todo proceso de análisis de datos, es necesario realizar una depuración y ajuste de los datos. Para este caso, fue necesario trabajar sobre la variable de la fecha y el valor.

# Vamos a sacar un valor promedio de cada variable por mes
#Para el dolar vamos a tomar el campo vigenciadesde

dolar <- dolar %>% 
  mutate(fecha = as.Date(vigenciadesde, format = "%d/%m/%y")) %>% 
  mutate(mes = month(fecha), anio = year(fecha)) %>% 
  mutate(valor = as.double(str_replace(valor,",",""))) %>% 
  group_by(anio,mes) %>%
  summarise(precio_dolar = mean(valor), .groups = "drop") 

A continuación se muestra una fracción de los datos depurados.

dolar %>%
  datatable()

Aguacate

hass %>% 
  head(5) %>% 
  ftable()

mercado

producto

fecha

precio_kg

Bogotá, D.C., Corabastos

Aguacate hass (Semanal | Mensual)

Sabado, Junio 23 de 2012

3,483

Bogotá, D.C., Corabastos

Aguacate hass (Semanal | Mensual)

Sabado, Junio 30 de 2012

3,459

Bogotá, D.C., Corabastos

Aguacate hass (Semanal | Mensual)

Sabado, Julio 7 de 2012

3,478

Bogotá, D.C., Corabastos

Aguacate hass (Semanal | Mensual)

Sabado, Julio 14 de 2012

3,493

Bogotá, D.C., Corabastos

Aguacate hass (Semanal | Mensual)

Sabado, Julio 21 de 2012

3,470

# 🖇️ Se genera el vector de meses para usarlo más abajo con la función match.
meses <- c("Enero", "Febrero", "Marzo", "Abril", "Mayo",
           "Junio", "Julio", "Agosto", "Septiembre", "Octubre",
           "Noviembre", "Diciembre")

hass <- hass %>% 
  # Otra forma de sacar el anio.
  #mutate(anio = substring(fecha,nchar(fecha)-4,nchar(fecha))) %>% 
  #mutate(espacio = grepl(", ",fecha))
  mutate(mes = sapply(strsplit(fecha, " "), "[", 2),
         anio = sapply(strsplit(fecha, " "), "[", 5)) %>% 
  mutate(anio = as.double(anio)) %>% 
  # 🖇 Se utiliza la función match con el vector meses.
  mutate(mes = match(mes,meses)) %>% 
  group_by(anio,mes) %>%
  summarise(precio_aguacate_kg = mean(precio_kg), .groups = "drop") 
hass %>% datatable()

Datos unidos

hass_dolar <- dolar %>% 
  right_join(hass, by = c("mes","anio")) 

hass_dolar %>% 
  head(5) %>% 
  ftable()

anio

mes

precio_dolar

precio_aguacate_kg

2,012

6

1,792.232

3,471.00

2,012

7

1,785.141

3,459.50

2,012

8

1,806.341

3,543.50

2,012

9

1,801.948

3,792.00

2,012

10

1,805.674

3,661.25

1.4 Diagrama de dispersión.

hass_dolar %>% 
  ggplot(aes(x= precio_dolar, y= precio_aguacate_kg )) +
  geom_point()+ theme_light()

Cómo se puede notar uno intuye que hay alguna relación entre ambas variables en la medida en que parece suceder que cuando el precio del dólar aumenta, también lo hace el precio de aguacate. La relación no es perfecta pues no vemos una curva suave que sea capaz de pasar por todos los puntos, más bien debemos reconocer que hay ciertas variaciones en el proceso. Estas variaciones pueden deberse a factores no contemplados por el modelo, después de todo no podemos esperar que el precio del aguacate dependa exclusivamente del precio del dólar.

1.5 Generando un modelo súper simplificado (modelo ingenuo)

Ahora vamos a generar el modelo más sencillo que podamos imaginar para predecir el precio del aguacate, al cual llamaremos un modelo ingenuo. La razón del apelativo es que el modelo “ingenuamente” supondrá que es posible predecir el precio que tendrá el aguacate con el promedio histórico.

La siguiente gráfica muestra la linea horizontal que simboliza el promedio histórico del precio del aguacate. También debe notar que predecir basado en esta línea es despreciar cualquier aporte de la variable precio del dólar en la predicción, es decir, es un modelo que no explota la relación (sea esta débil o sea esta fuerte) entre el precio del aguacate y su variable explicativa precio del dólar.

ggplot(data = hass_dolar,
       aes(x= precio_dolar, y= precio_aguacate_kg )) +
geom_point()+ 
geom_hline(yintercept = mean(hass_dolar$precio_aguacate_kg),
           color = "red")+
geom_text(data=hass_dolar,
          aes(x= 5700, y = mean(hass_dolar$precio_aguacate_kg),
              label = paste("Precio promedio histórico: $",
                            round(mean(hass_dolar$precio_aguacate_kg), 2))),
          hjust = 1.2, vjust = -0.2, color = "red"
) +
theme_light()

Observe que el modelo ingenuo no sería tan inadecuado en dos circunstancias:

  • Si los datos se encuentran muy cercanos al precio histórico de manera consistente sin importar en que valor del precio del dólar me posicione. Lo que vendría a indicar un escenario donde moverme a lo largo de diferentes valores del precio del dólar no genera ningún patrón de distanciamiento frente al valor histórico de referencia.

  • Si los datos del precio del aguacate se alejan de la línea histórica de referencia pero el precio del dólar tampoco puede seguirlos, es decir, si la relación entre precio del dólar y precio del aguacate se parece a una nube de puntos sin ningún patrón o tendencia evidente. En este caso cualquier cosa que supongamos de la relación entre precio de dólar y precio del aguacate será producto de nuestra imaginación 🤦🤦

En estos dos escenarios es mejor retener el modelo ingenuo a falta de una variable que de verdad explique lo que pasa con el precio del aguacate.

Ahora bien, ninguno de los dos escenarios anteriores parece ser el nuestro, más bien acá se nota que el precio del dólar si está acompasado con el precio del aguacate. Sin embargo antes de entrar en la búsqueda de esas relaciones es importante notar dos cosas:

  • El modelo ingenuo es mi modelo de referencia, esto significa que si un modelo predictivo construido para predecir el precio del aguacate es mejor, lo tendrá que ser respecto a este modelo ingenuo.

  • Es posible medir el desajuste actual de mi modelo ingenuo tomando la distancia de cada punto a la recta horizontal de promedio histórico, esto es:

1.6 Calidad de ajuste del modelo ingenuo

Sea \(e_i = y_i - \bar{y}\) las diferencia entre el dato del precio de aguacate \(y_i\) y el promedio histórico \(\bar{y}\). Como tengo muchas diferencias, será necesario definir una métrica resumen, por conveniencia elijamos la suma de los cuadrados de las diferencias, debido a que si sumamos las diferencias, éstas se me compensarán, pues diferencias negativas cancelarán las positivas, y no quiero esto. Más bien me interesa que ambas sumen a mi medida de desajuste. Llamemos a esta medida la suma de cuadrados del error:

\[SCE_{\text{ingenuo}} = \sum^{n}_{i=1}e_i²=\sum^{n}_{i=1}(y_i - \bar{y})²\]

En la definición de la cantidad quisimos usar el subíndice “ingenuo” para nombrar al SCE, esto con el fin de hacer énfasis que tal métrica se asocia al modelo, si el modelo cambia a otro tipo de modelo el SCE cambiará también: será la distancia de cada \(y_i\) a la curva definida por el otro modelo. Veremos esto más adelante. No obstante ya podemos sacar nuestras primeras conclusiones:

  1. Entre más grandes sean cada una de las distancias \(e_i\) individuales, más grande será SCE.
  2. Todas las distancias contribuyen con igual importancia al SCE, excepto en lo que se refiere a su magnitud no hay porque pensar que un distanciamiento de \(u\) unidades sea más importante que otro de las mismas \(u\) unidades. Esto puede ser obvio pero no lo es cuando, por ejemplo, queremos castigar más las distancias que se dan en la zona central del gráfico que las que se dan en los extremos del gráfico.
  3. El SCE es sensible a la cantidad de datos, razón por la cual entre más datos (es decir entre más grande sea n), más grande será SCE. Esto dificulta seriamente la posibilidad de comparar modelos que fueron calculados sobre una cantidad distinta de puntos.

De acuerdo con lo anterior modifiquemos un poco nuestra definición de desajuste, así:

\[MSE_{\text{ingenuo}} = \color{red}{\frac{1}{n}}\sum^{n}_{i=1}(y_i - \bar{y})²\]

Es decir, nuestra medida de desajuste será el promedio de la suma de las diferencias cuadradas de cada dato \(y_i\) respecto a su media \(\bar{y}\), dicho de otra manera la varianza de la variable \(Y\) a predecir!!!!

A continuación invitamos al estudiante a calcular la varianza de la variable a predecir, es decir, el desajuste del modelo ingenuo.

# En esta sección el estudiante desarrollará los cálculos.

# 💡 Recordemos que la varianza de los residuales respecto del modelo ingenuo, es la misma varianza de la variable a predecir.

Observe el siguiente gráfico en el que hemos pintado las distancias que hay entre cada dato \(y_i\) observado y el promedio histórico (valor predicho \(\bar{y}\) por el modelo ingenuo)

hass_dolar %>% 
  ggplot(aes(x= precio_dolar, y= precio_aguacate_kg )) +
  geom_point()+ 
  geom_hline(yintercept = mean(hass_dolar$precio_aguacate_kg))+
  geom_segment(aes(xend=precio_dolar, yend=mean(precio_aguacate_kg)),
               col='red', lty='dashed')+
  theme_light()

1.7 Construyendo un modelo lineal simple con los datos

Este es el momento de comenzar a usar la información extra con la que contamos, esto es, el precio del dólar, el cual se convertirá en una variable predictora del precio del aguacate. En estadística no hay una única forma de hacer esto, pero lo usual es partir de especificaciones muy sencillas. Veamos:

Sea \(Y\) el precio del aguacate y sea \(X\) el precio del dólar. Investiguemos si este modelo sencillo nos funciona:

\[y = \beta_{0} + \beta_{1} * x + e \text{ con } e \text{~}N(0, \sigma^2)\]

1.8 Algunas observaciones sobre las fortalezas y limitaciones del modelo

Expliquemos qué implica la ecuación anterior haciendo ciertas anotaciones de gran relevancia práctica para entender cómo los estadísticos piensan al momento de plantear modelos:

Observación genial 1:

Para un x fijo, es decir un x dado, es posible calcular el valor de \(y\) con el modelo anterior, donde tácitamente estamos diciendo que el valor de \(y\) dependerá de \(x\) en forma exacta, excepto por una perturbación aleatoria \(e\) que representará la parte del valor de \(y\) que no puede ser capturada por el término \(\beta_{0} + \beta_{1} * x\).

¿Cuándo sucederá que el valor de \(y\) pueda ser exactamente capturado por el término \(\beta_{0} + \beta_{1} * x\)? Cuando \(y\) dependa en forma lineal exacta del valor de \(x\). Esta es una condición demasiado fuerte que rara vez ocurre en la práctica, por eso agregamos la perturbación \(e\) como una forma de modelar la incertidumbre en la determinación de \(y\) dado un valor de \(x\). Esto no soluciona todos los problemas pero ayuda a obtener un modelo relativamente más realista.

¿Cuáles son los riesgos que implica asumir a priori un comportamiento específico para el error?

Observación genial 2:

La perturbación \(e\) puede recibir varios nombres en estadística, la encontrarás nombrada como: perturbación, choque, error, o residual. Lo importante no es el nombre que reciba, sino las condiciones que vamos a suponer sobre sus valores.

En particular no vamos a exigir que \(e\) tenga valores predecibles, muy al contrario vamos a exigir que sea un valor aleatorio, precisamente porque está modelando incertidumbre. Pero esto no implica que no podamos poner ciertas condiciones sobre el tipo de aleatoriedad deseada para \(e\). En este caso supondremos que \(e\) se distribuye como una variable aleatoria normal con media \(\mu=0\) y desviación estándar \(\sigma\). Conviene aclarar que esto es una suposición, no implica que efectivamente los errores vayan a cumplir tal supuesto.

Observación genial 3:

¿Por qué suponemos normalidad para \(e\)? Porque si \(e\) es en efecto una especie de componente incierto en el valor de \(y\), entonces muy seguramente será la suma de muchos efectos independientes que lo están provocando. En nuestro contexto estos efectos inciertos pueden ser:

  • \(F_1\): Cambios en la productividad de los cultivos.
  • \(F_2\): Condiciones climatológicas.
  • \(F_3\): Precios de los insumos agrícolas.
  • \(F_4\): Capacidad de negociación de los productores.
  • \(F_5\): Presencia o no de subsidios estatales.
  • \(F_6\): Precios de productos complementarios o sustitutos.

  • \(F_k\): precio de la gasolina.

Y un gran etcétera.

Es decir, existen innumerables factores, distintos al precio del dólar que impiden que el precio del aguacate pueda ser determinado de forma exacta sólo por observar el valor del dólar. En estadística se ha estudiado que cuando no estamos midiendo los demás factores que afectan el valor de una variable, y dejamos que estos contribuyan de forma desacoplada e independiente a dicho valor, el efecto general \(e\) que resume el efecto global de todos los factores a la vez, suele comportarse bajo la distribución normal sencillamente porque hay un teorema en matemáticas, conocido como el Teorema del Límite Central que básicamente así lo garantiza, al postular que: la suma de \(k\) efectos aleatorios independientes tiende a adoptar el comportamiento normal conforme la cantidad k de efectos crece.

Es importante destacar que el Teorema del Límite Central tiene algunas condiciones y suposiciones, como la independencia de las variables aleatorias y que la suma debe efectuarse sobre una cantidad \(k\) de ellas suficientemente grande. Éste es uno de los conceptos fundamentales en estadística y es ampliamente utilizado en la teoría y la práctica de esta disciplina.

Observación para nada genial 4 😥😢:

Excepto por el componente aleatorio \(e\), el modelo restringe la relación de \(Y\) y \(X\) al universo de relaciones lineales. Existirá una posible relación por cada par de valores \(\beta_0\) y \(\beta_1\) que decidamos elegir. Pero por más que nos esforcemos en modificarlos siempre conducirán a relaciones representadas por líneas rectas, de ahí que el modelo asuma el nombre de regresión lineal.

Si queremos capturar otras posibles dependencias de \(Y\) respecto a \(X\) podemos generalizar la relación usando \(y = f(x) + e\) donde \(f(x)\) puede ser una nueva especificación funcional con otra estructura deseada, por ejemplo un polinomio o cualquier otra función que querramos definir. Lo importante es que detrás de la definición de \(f(x)\) haya alguna justificación producto de haber analizado los datos en búsqueda de relaciones que capturen bien la dependencia.

Observación genial final:

Son los datos observados para cada valor de \(y\) y \(x\) los que nos deben informar sobre el par de valores \(\beta_0\) y \(\beta_1\) que crean la relación lineal que mejor representa nuestra nube de puntos, pero de nada sirven los valores observados si no definimos una métrica que nos informe sobre la calidad del ajuste.

1.9 Introducción a las medidas de ajuste para modelos no ingenuos

De la observación final del apartado anterior nos queda una inquietud: ¿Cómo decidimos el par de valores a asignar a los parámetros del modelo lineal si no tenemos un criterio para poder saber cuál es el mejor modelo?

Para salir de este embrollo, reconozcamos que nunca en un problema real sujeto a incertidumbre una recta pasará exactamente por todos los puntos, razón por la cual aparecerán componentes de error \(e_i\) por cada \(y_i\) y \(x_i\) observado. Vimos que esto ocurrión con el modelo ingenuo y por supuesto ocurre también para un modelo no ingenuo. La idea es entonces reducir la magnitud de estos errores al mínimo posible y el par de valores \(\beta_1\) y \(\beta_2\) que así lo logren será nuestra elección óptima de cara al objetivo de minimización del error.

Otro criterio de ajuste que se podría usar es el de maximizar la probabilidad de ocurrencia de los valores observados bajo el supuesto de que el modelo usa valores \(\beta_0\) y \(\beta_1\) previamente especificados. De esta manera si un par de valores hace menos probable haber obtenido nuestros datos observados deberán descartarse.

¿Pero cómo podemos medir la probabilidad de ocurrencia de nuestros datos observados una vez damos valores específicos para los betas? Esta es una cuestión que se tratará más adelante cuando discutamos los métodos de estimación de parámetros, nombre con el que se conoce al procedimiento estadístico encaminado a elegir los mejores betas para un conjunto de datos observados.

Por lo pronto vamos a proceder de manera más intuitiva, dejando de lado el formalismo matemático, y vamos a ejemplificar un posible modelo ajustado usando valores \(\beta_0=700\) y \(\beta_1 = 1.2\). Hemos elegido estos valores por simple inspección visual así que no esperamos que sean óptimos en ningún sentido, veamos:

# Calculamos las predicciones de "y" usando el modelo
b0 = 700
b1 = 1.2
x= hass_dolar$precio_dolar
y_pred = b0 + b1*x

hass_dolar %>% 
  ggplot(aes(x= precio_dolar, y= precio_aguacate_kg ))+
  geom_point() +
  geom_line( aes( x= x, y = y_pred), col = "red")

Aunque no estamos seguros de la optimalidad de los betas elegidos, el ajuste basado en este modelo es a simple vista decente, pero no deberíamos depender de esto. Es necesario usar las herramientas que R nos ofrece para estimar los betas. Antes de pasar a su uso ciego, revisemos en qué consideraciones estadísticas se basan estos métodos. Para ello definamos primero los errores cometidos por el modelo, y más aún, la medida de desempeño que usaremos para decidir cuál es la mejor elección de los betas. Usaremos entonces la métrica MSE(Mean Squared Errors) dada por:

\[MSE_{\text{NO ingenuo}} = \frac{1}{n}\sum^{n}_{i=1}e_i²\]

Por otro lado podemos decir que \(e_i = y_i - \hat{y}\), donde la cantidad \(\hat{y}\) representa el valor de la variable \(Y\) que arroja el modelo, el cual es por regla general diferente al verdadero valor observado de \(Y\), es decir, los errores son las diferencia de lo observado \(y_i\) y lo predicho por el modelo basado en la recta \(\hat{y}\), por lo tanto:

\[MSE_{\text{NO ingenuo}} = \frac{1}{n}\sum^{n}_{i=1}(y_i - \hat{y})²\]

Y más específicamente:

\[MSE_{\text{NO ingenuo}} = \frac{1}{n}\sum^{n}_{i=1}(y_i - (\beta_0 + \beta_1*x_i))²\] \[MSE_{\text{NO ingenuo}} = \frac{1}{n}\sum^{n}_{i=1}(y_i - \beta_0 - \beta_1*x_i)²\]

Con esta expresión ya se hace patente que \(MSE_{\text{NO ingenuo}}\) es función de los pares de valores \((x_i, y_i)\) a los que se debe ajustar la recta, así como de los valores que elijamos para los betas.

Sin embargo, como los pares \((x_i, y_i)\) ya vienen dados por nuestra base de datos o, hablando en términos más generales, por la información recaudada para la construcción del modelo, sólo tendremos opción de variar los betas en aras de minimizar la cantidad \(MSE_{\text{NO ingenuo}}\). Dicho de otra manera los pares de valores \((x_i, y_i)\) se comportan como constantes a lo largo de mi proceso de minimización.

Dado lo anterior, es buena idea definir una función MSE en R que permita calcular, para diferentes elecciones de los betas, la medida de desempeño. Esta medida de desempeño es conocida también en la literatura como función de pérdida, queriendo indicar con ello una pérdida de ajuste del modelo a los datos. Es importante aclarar que no existe una única elección para la función de pérdida de un modelo predictivo, otras alternativas típicas son la mediana de las desviaciones absolutas (MEDA), la media de los errores absolutos (MAE), entre otros. Sin embargo en regresión lineal es usual usar MSE porque su expresión como cuadrado de errores permite usar teoría de inferencia basada en la distribución F cuando los errores se asumen normales. Veremos esto más a detalle cuando discutamos pruebas de hipótesis sobre el modelo de regresión.

1.10 Una relación importante entre los betas del modelo lineal

Retomando nuestro modelo original

\[y = \beta_{0} + \beta_{1} * x + e \text{ con } e \text{~}N(0, \sigma^2)\]

Tomando valor esperado en ambos lados de la ecuación tenemos que:

\[E(y) = E(\beta_{0} + \beta_{1} * x + e)\] Pero como \(\beta_0\) y \(\beta_1\) son constantes, tenemos que:

\[E(y) = E(\beta_0) + E(\beta_1*x) + E(e)\] Así que:

\[E(y) = \beta_0 +\int_{-\infty}^{\infty}\beta_1 *xf(x)dx + E(e)\] \[E(y) = \beta_0 + \beta_1*\int_{-\infty}^{\infty}xf(x)dx + E(e)\] Y como \(\int_{-\infty}^{\infty}xf(x)dx\) es por definición \(E(x)\) y además \(E(e) = 0\) por diseño (de acuerdo con el supuesto usado para la distribución del error), entonces:

\[E(y) = \beta_0 + \beta_1*E(x)\] Dicho de otra manera, la recta del modelo debe pasar por el punto \((E(x), E(y))\) !!

La recta debe pasar por el centroide de la nube de puntos.

Ciertamente no conocemos el valor esperado ni de la variable predictora X ni de la variable respuesta Y, pero buenos estimados basados en información muestral serían \(\bar{x}\) y \(\bar{y}\). Es decir, no podremos saber con precisión los valores \(\beta_0\) y \(\beta_1\), pero nos podemos conformar con buenas estimaciones que cumplan la condición:

\[\bar{y} = \color{red}{\widetilde{\beta}_0} + \color{red}{\widetilde{\beta}_1}\bar{x}\]

El coloreado y la tilde ancha para los betas pretenden enfatizar que estos valores son estimados de los correspondientes \(\beta_0\) Y \(\beta_1\) del modelo teórico subyacente. En conclusión: la mejor recta estimada deberá pasar por \((\bar{x}, \bar{y})\), y por lo tanto hay una relación atando a los betas estimados:

\[\begin{equation} \label{eq:ec0} \color{red}{\widetilde{\beta}_0} = \bar{y} - \color{red}{\widetilde{\beta}_1}\bar{x} \end{equation} \] Es decir que dando valores a \(\color{red}{\widetilde{\beta}_1}\) podremos obtener valores para \(\color{red}{\widetilde{\beta}_0}\) y en últimas dibujar diversas rectas en el plano, con el fin de establecer cuál de ellas minimiza el \(MSE_{\text{NO ingenuo}}\).

1.11 Una animación que ilustra la estimación de betas

A continuación presentamos un código en R que permite comprender cómo funcionan numéricamente los métodos de estimación de betas para un modelo lineal simple, específicamente cuando se usa como método de estimación el criterio de minimización del promedio de los cuadrados de los errores:

# Extraemos las variables de interés
x = hass_dolar$precio_dolar
y = hass_dolar$precio_aguacate_kg

# Calculamos medias muestrales para las variables X y Y

x_barra <- mean(x)
y_barra <- mean(y)

# Definimos posibles valores para b1, y por consiguiente para b0
b1_values <- seq(-1, 3, by = 0.08)
b0_values <- y_barra - b1_values * x_barra

# Fijamos límites en X y Y para cada gráfico
xmin = min(x) - 50
xmax = max(x) + 50
ymin = min(y) - 50
ymax = max(y) + 50

# Construimos una función para calcular el MSE en función de x_i, y_i, b1 y b0
MSE <- function(x, y, b0, b1) {
  errors <- y - (b0 + b1*x)
  squared_errors <- errors^2
  mse <- mean(squared_errors)
  return(mse)
}

# Crear un dataframe para almacenar los resultados
list_mse <- vector()

# Calcular el MSE para diferentes valores de B1
for (i in 1:length(b1_values)) {
  mse <- MSE(x, y, b0_values[i], b1_values[i])
  list_mse[i] <- mse
}

mse_df = data.frame(
  B0 = b0_values,
  B1 = b1_values,
  mse = list_mse,
  color_point = 'green'
)

# Capturando el renglón donde ocurre el mínimo mse
pos_min <- which(mse_df$mse==min(mse_df$mse))

# Impriendo las filas cercanas al mínimo
mse_df$color_point[(pos_min - 4):(pos_min + 4)] <- rep('red',9)
ftable(mse_df[(pos_min - 8):(pos_min + 8), ])

B0

B1

mse

color_point

3,475.69602

0.28

584,101.9

green

3,230.39243

0.36

516,278.7

green

2,985.08884

0.44

457,137.2

green

2,739.78525

0.52

406,677.3

green

2,494.48166

0.60

364,899.1

red

2,249.17807

0.68

331,802.7

red

2,003.87448

0.76

307,387.9

red

1,758.57089

0.84

291,654.8

red

1,513.26730

0.92

284,603.3

red

1,267.96371

1.00

286,233.6

red

1,022.66012

1.08

296,545.5

red

777.35653

1.16

315,539.2

red

532.05294

1.24

343,214.5

red

286.74934

1.32

379,571.5

green

41.44575

1.40

424,610.2

green

-203.85784

1.48

478,330.6

green

-449.16143

1.56

540,732.6

green

¿Qué hemos logrado hasta acá?

  • Demostrar que el valor de \(MSE_{\text{NO ingenuo}}\) depende de los betas elegidos para el modelo.
  • Demostrar cómo el valor de \(\color{red}{\widetilde{\beta}_0}\) está atado con el de \(\color{red}{\widetilde{\beta}_1}\) si de verdad queremos que la recta estimada respete la condición de que \(e \text{~}N(0, \sigma^2)\) y por lo tanto de que \(E(e)=0\).
  • Generar un dataframe que calcula los valores de \(MSE_{\text{NO ingenuo}}\) para cada par de valores elegidos para \(\color{red}{\widetilde{\beta}_0}\) y \(\color{red}{\widetilde{\beta}_1}\), dando valores a \(\color{red}{\widetilde{\beta}_1}\)
  • Estimamos que el \(\color{red}{\widetilde{\beta}_1}\) óptimo está cerca del valor 0.94, y que el \(\color{red}{\widetilde{\beta}_0}\) óptimo está cerca del valor 1452.

Ahora graficaremos una curva que muestre esta relación:

if (!file.exists("beta_vs_mse.gif")){
  # Crear la animación
  animation <- ggplot(mse_df, aes(B1, mse)) +
    geom_line(aes( x= B1, y = mse, color = color_point), linetype = "dashed") +
    geom_point(aes(color = color_point), size = 3) +
    labs(title = "Valores de MSE para cada posible B1, respetando E(e) = 0") +
    geom_text(data = mse_df, 
              aes(label = paste("MSE: ", round(mse,0))), x= 1, y = 2e6, color='black'
    ) +
    geom_text(data = mse_df, 
              aes(label = paste("MSE mínimo: ", round(min(mse),0))),
              x= 1, y = min(mse_df$mse) - 1e5, color='red') +
    scale_color_manual(values = c("red" = "red", "green" = "green"),
                       labels = c("Subóptima", "Óptima"),
                       name = "Soluciones") +
    theme_minimal() +
    transition_reveal(B1)

  animate(animation, duration = 10, fps = 2, renderer = gifski_renderer())
  anim_save("beta_vs_mse.gif")
}

Img 1: Relación entre B1 y MSE

Ahora veamos el efecto que tiene la elección de los betas en la recta regresora, e incluyamos pausas en la animación para detenernos en el momento en que encontremos la recta óptima.

if (!file.exists("beta_estimation.gif")){
  
# Las pausas en la animación son fotogramas con información repetida correspondiente
# al renglón óptimo

fotogramas_pausa <- data.frame(
  B0 = rep(b0_values[pos_min], 2*nrow(mse_df)),
  B1 = rep(b1_values[pos_min], 2*nrow(mse_df)),
  mse = rep(list_mse[pos_min], 2*nrow(mse_df)),
  color_point = rep('red',2*nrow(mse_df))  
)

mse_df <- rbind(mse_df, fotogramas_pausa)
mse_df <- mse_df[order(mse_df$B1), ]
mse_df$estado <- seq(1,nrow(mse_df))

# Construimos el dataframe de predicciones para cada estado, con el cual poder
# graficar los residuales en cada fotograma

df_predicciones <- data.frame(B0_pred=numeric(0),
                              B1_pred=numeric(0),
                              x_from_pred=numeric(0),
                              y_from_pred=numeric(0),
                              y_pred_mod=numeric(0),
                              estado = numeric(0))
for (i in 1:nrow(mse_df)){
  df_nuevo <-data.frame(
    BO_pred = rep(mse_df[i, 'B0'], length(x)),
    B1_pred = rep(mse_df[i, 'B1'], length(x)),
    x_from_pred = x,
    y_from_pred = y,
    y_pred_mod = mse_df[i, 'B0'] + mse_df[i,'B1']*x,
    estado = rep(mse_df[i, 'estado'], length(x))
  ) 
  df_predicciones <- rbind(df_predicciones, df_nuevo)
}

# Creamos el gráfico base

base_plot <- ggplot(data.frame(x = x, y = y), aes(x = x, y = y)) +
  geom_point() +
  geom_point(x=mean(x), y=mean(y), color='green', size=2.5) +
  geom_abline(aes(slope = 0, intercept = y_barra), 
              linetype = "dashed", color = "blue") +
  labs(title = "Estimación de Betas por Minimización del MSE") +
  theme_minimal()

# Creamos la animación

animation <- base_plot +
  geom_abline(aes(slope = B1, intercept = B0, color = color_point), data=mse_df) +
  geom_text(aes(label = paste("MSE: ", round(mse,0))), x= 2000, y = 6000, data=mse_df,
            color = "red")+
  geom_segment(data = df_predicciones,
               aes(x=x_from_pred, y=y_from_pred, xend=x_from_pred, yend=y_pred_mod),
               col = "violet", lty='dashed') +
  transition_states(estado, transition_length = 2, state_length = 1) +
  scale_color_manual(values = c("red" = "red", "green" = "green"),
                     labels = c("Subóptima", "Óptima"),
                     name = "Rectas regresoras") +
  enter_fade() +
  exit_fade()

# Guardamos la animación en un archivo GIF
animate(animation, duration = 20, fps = 2, renderer = gifski_renderer())
anim_save("beta_estimation.gif", animation)
}

Img 1: Estimación de betas por minimización de MSE

1.12 Usando funciones en R para estimar parámetros

Finalmente hemos llegado al punto en que usaremos las funciones de R para estimar parámetros del modelo de regresión lineal. Esto se hará bajo el enfoque de optimización, para ello conviene darle una mirada al siguiente problema:

Estime un par de valores \((\widetilde{\beta}_0,\widetilde{\beta}_1)\) tal que minimice la expresión \(MSE = f(\beta_0,\beta_1)\), donde:

\[f(\beta_0,\beta_1) = \frac{1}{n}\sum^{n}_{i=1}(y_i - \beta_0 - \beta_1*x_i)²\] Tomando en cuenta que no hay restricciones aplicadas a los betas.

Para resolver el problema primero efectuemos algunas simplificaciones tomando partida de la siguiente identidad algebraica: \((a+b+c)² = a²+b²+c²+2ab+2ac+2bc\), de tal manera que:

\[f(\beta_0,\beta_1) = \frac{1}{n}\sum^{n}_{i=1}(y_i²+\beta_0²+\beta_1²x_i²-2y_i\beta_0-2y_i\beta_1x_i+2\beta_0\beta_1x_i)\]

Ahora distribuyamos la suma y separemos aparte los términos que dependen de su índice:

\[f(\beta_0,\beta_1) = \frac{1}{n}\left[\sum^{n}_{i=1}y_i²+\beta_0²\sum^{n}_{i=1}1+\beta_1²\sum^{n}_{i=1}x_i²-2\beta_0\sum^{n}_{i=1}y_i-2\beta_1\sum^{n}_{i=1}y_ix_i+2\beta_0\beta_1\sum^{n}_{i=1}x_i\right]\] Simplificando términos y distribuyendo la fracción:

\[f(\beta_0,\beta_1) = \frac{1}{n}\sum^{n}_{i=1}y_i²+\frac{\beta_0²}{n}n+\frac{\beta_1²}{n}\sum^{n}_{i=1}x_i²-\frac{2\beta_0}{n}\sum^{n}_{i=1}y_i-\frac{2\beta_1}{n}\sum^{n}_{i=1}y_ix_i+\frac{2\beta_0\beta_1}{n}\sum^{n}_{i=1}x_i\] Para finalmente obtener:

\[f(\beta_0,\beta_1) = \beta_0² +\bar{x²}\beta_1² -2\bar{y}\beta_0-\left[\frac{2}{n}\sum^{n}_{i=1}x_iy_i\right]\beta_1+2\bar{x}\beta_0\beta_1+\bar{y²} \] En la expresión anterior, las cantidades \(\bar{x²}\), \(-2\bar{y}\), \(-\frac{2}{n}\sum^{n}_{i=1}x_iy_i\), \(2\bar{x}\), \(\bar{y²}\) son todas constantes por depender sólo de los datos que alimentan al modelo, en ese sentido son coeficientes numéricos para los respectivos betas de la función.

LLegados a este punto, podemos sacar las siguientes conclusiones:

  • La función \(f\) sólo depende de los betas, puesto que los demás valores son constantes, es decir los valores \((x_i, y_i)\) son números consignados en una base de datos, y \(n\) es la cantidad de pares de números consignados. Todos los coeficientes dependerán de operaciones sobre estos números.

  • De acuerdo con lo anterior la función \(f\) es una función multivariada, pues depende de dos variables: \(\beta_0\) y \(\beta_1\).

  • Existen métodos muy sencillos para resolver un problema de optimización no restringido cuando la función a optimizar es convexa, entonces sólo nos resta demostrar que la función anterior lo es.

Definición de convexidad

Para demostrar que \(f(\beta_0, \beta_1)\) es convexa, debemos verificar que su matriz hessiana sea semidefinida positiva. La matriz hessiana es una matriz de segundas derivadas parciales de la función.

Primero, calculemos las derivadas parciales de \(f\) con respecto a \(\beta_0\) y \(\beta_1\):

\[\begin{equation} \label{eq:ec1} \frac{\partial f}{\partial \beta_0} = 2\beta_0-2\bar{y}+2\bar{x}\beta_1 \end{equation} \] \[\begin{equation} \label{eq:ec2} \frac{\partial f}{\partial \beta_1} = 2\bar{x²}\beta_1-\frac{2}{n}\sum^{n}_{i=1}x_iy_i+2\bar{x}\beta_0 \end{equation} \] Luego, calculemos las segundas derivadas parciales de \(f\) con respecto a \(\beta_0\) y \(\beta_1\):

\[\frac{\partial^2 f}{\partial \beta_0^2} = 2 > 0\] \[\frac{\partial^2 f}{\partial \beta_1^2} = 2\bar{x^2} > 0\] \[\frac{\partial^2 f}{\partial \beta_0 \partial \beta_1} = 2\bar{x}\]

La matriz hessiana de \(f\) es:

\[ H(f) = \begin{bmatrix} \frac{\partial^2 f}{\partial \beta_0^2} & \frac{\partial^2 f}{\partial \beta_0 \partial \beta_1} \\ \frac{\partial^2 f}{\partial \beta_0 \partial \beta_1} & \frac{\partial^2 f}{\partial \beta_1^2} \end{bmatrix} \]

Por lo tanto:

\[ H(f) = \begin{bmatrix} 2 & 2\bar{x}\\ 2\bar{x} & 2\bar{x^2} \end{bmatrix} \]

Para que \(f(\beta_0, \beta_1)\) sea convexa, la matriz hessiana debe ser semidefinida positiva, lo que significa que debe cumplir la propiedad: \[v^ \top Hv>=0 \text{ para }\forall v \in \mathbb{R}² \] Es decir, debemos verificar si:

\[\begin{bmatrix} a & b \end{bmatrix}H\begin{bmatrix} a \\ b \\ \end{bmatrix}>=0 \text{ para } \forall (a,b) \]

Veamos si esto se cumple:

\[ \begin{bmatrix} a & b \end{bmatrix} \begin{bmatrix} 2 & 2\bar{x} \\ 2\bar{x} & 2\bar{x^2} \end{bmatrix} \begin{bmatrix} a \\ b \\ \end{bmatrix}\\= \begin{bmatrix} 2a+2b\bar{x} & 2a\bar{x} + 2b\bar{x²} \end{bmatrix} \begin{bmatrix} a \\ b \\ \end{bmatrix}\\=2a²+2ab\bar{x}+2ab\bar{x}+2b²\bar{x²} \\=2a²+4ab\bar{x}+2b²\bar{x²}\\ =2(a²+2ab\bar{x}+b²(\bar{x})²)-2b^2(\bar{x})²+2b²\bar{x²}\\ =2(a+b\bar{x})²+2b²(\bar{x²}-(\bar{x})²) \]

Acá es interesante notar que:

\[\bar{x²}-(\bar{x})²=\frac{1}{n}\left[\sum^{n}_{i=1}x_i²\right]-\bar{x}\bar{x} =\frac{1}{n}\left[\sum^{n}_{i=1}x_i²-\bar{x}(n\bar{x})\right]\\ =\frac{1}{n}\left[\sum^{n}_{i=1}x_i²-\bar{x}\sum^{n}_{i=1}x_i\right] =\frac{1}{n}\left[\sum^{n}_{i=1}(x_i²-\bar{x}x_i)\right] \\ =\frac{1}{n}\left[\sum^{n}_{i=1}((x_i²-2\bar{x}x_i)+\bar{x}x_i)\right] =\frac{1}{n}\left[\sum^{n}_{i=1}\left((x_i²-2\bar{x}x_i+\bar{x}²)+(\bar{x}x_i-\bar{x}²)\right)\right]\\ =\frac{1}{n}\left[\sum^{n}_{i=1}\left((x_i-\bar{x})²+\bar{x}(x_i-\bar{x})\right)\right] =\frac{1}{n}\sum^{n}_{i=1}(x_i-\bar{x})²+\bar{x}\left[\sum^{n}_{i=1}\frac{x_i}{n}-\sum^{n}_{i=1}\frac{\bar{x}}{n}\right]\\ =\frac{S_{xx}}{n}+\bar{x}\left(\bar{x}-\bar{x}\sum^{n}_{i=1}\frac{1}{n}\right)\\ =\frac{S_{xx}}{n}+\bar{x}(\bar{x}-\bar{x}*1)=\frac{S_{xx}}{n}\]

Obteniendo finalmente

\[\begin{equation} \label{eq:sigma_x} \bar{x²}-(\bar{x})² = \frac{S_{xx}}{n} \end{equation} \] Note que aquí hemos usado la siguiente definición:

\[S_{xx}=\sum^{n}_{i=1}(x_i-\bar{x})*(x_i-\bar{x})=\sum^{n}_{i=1}(x_i-\bar{x})^2\] Luego:

\[v^ \top Hv = 2(a+b\bar{x})²+2b²\frac{S_{xx}}{n} \text{ con } \text{ }v=(a,b)\] Y como ambos sumandos son no negativos se tiene que \[v^ \top Hv >= 0 \text{ para }\forall v\in \mathbb{R}²\] Es decir, la matriz hessiana es semidefinida positiva y por lo tanto \(f(\beta_0,\beta_1)\) es convexa en el espacio de parámetros \(\beta_0\) y \(\beta_1\). Esto implica que el problema de regresión lineal, que busca minimizar esta función de errores, tiene un mínimo global y puede ser resuelto de manera eficiente mediante métodos de optimización convexa.

Una condición necesaria y suficiente para que \((\widetilde{\beta}_0,\widetilde{\beta}_1)\) sea un mínimo global es que \(\nabla{} f(\widetilde{\beta}_0,\widetilde{\beta}_1)=0\)

Por lo tanto retomando \(\ref{eq:ec1}\) y \(\ref{eq:ec2}\) el problema de minimización queda resuelto hallando la solución al sistema de ecuaciones:

\[\begin{equation} \label{eq:ec3} 2\widetilde{\beta}_0-2\bar{y}+2\bar{x}\widetilde{\beta}_1=0 \end{equation} \] \[\begin{equation} \label{eq:ec4} 2\bar{x²}\widetilde{\beta}_1-\frac{2}{n}\sum^{n}_{i=1}x_iy_i+2\bar{x}\widetilde{\beta}_0=0 \end{equation} \]

De modo que despejando \(\widetilde{\beta}_0\) en \(\ref{eq:ec3}\) tenemos lo siguiente:

\[\begin{equation} \label{eq:ec5} \widetilde{\beta}_0 = \bar{y} - \widetilde{\beta}_1\bar{x} \end{equation} \] Esto en concordancia con \(\ref{eq:ec0}\), y además haciendo reemplazos en \(\ref{eq:ec4}\) obtenemos:

\[ 2\bar{x²}\widetilde{\beta}_1-\frac{2}{n}\sum^{n}_{i=1}x_iy_i+2\bar{x}(\bar{y}-\widetilde{\beta}_1\bar{x})=0\\ 2\bar{x²}\widetilde{\beta}_1-2\bar{x}²\widetilde{\beta}_1 = \frac{2}{n}\sum^{n}_{i=1}x_iy_i-2\bar{x}\bar{y}\\ (\bar{x²} - \bar{x}²)\widetilde{\beta}_1 = \frac{1}{n}\sum^{n}_{i=1}x_iy_i-\bar{x}\bar{y}\\\]

Y usando la expresión dada en \(\ref{eq:sigma_x}\) tenemos:

\[\frac{S_{xx}}{n}\widetilde{\beta}_1=\frac{1}{n}\left(\sum^{n}_{i=1}x_iy_i-\frac{\bar{x}\bar{y}}{n}\right)\\S_{xx}\widetilde{\beta}_1 = \sum^{n}_{i=1}x_iy_i-\frac{\bar{x}\bar{y}}{n} \]

Ahora, usando el hecho de que el lado derecho de la ecuación anterior cumple la siguiente equivalencia:

\[\sum^{n}_{i=1}x_iy_i-\frac{\bar{x}\bar{y}}{n}= \sum^{n}_{i=1}(x_i-\bar{x})*(y_i-\bar{y}) = S_{xy}\] Es decir, se convierte en la suma de productos cruzados de \(x\) y la variable respuesta \(y\), una vez estos han sido centrados respecto a la media, así como :

\[\begin{equation} \label{eq:ec6} \widetilde{\beta}_1 = \frac{S_{xy}}{S_{xx}} \end{equation} \]

Las ecuaciones \(\ref{eq:ec5}\) y \(\ref{eq:ec6}\) son muy informativas en su estructura, la primera indica que la recta regresora debe pasar por el centroide de la nube de puntos: \((\bar{x}, \bar{y})\) y la segunda indica que la pendiente de la recta es proporcional al grado de relación lineal que tenga la variable \(x\) con la variable \(y\) medida a través de la suma de productos cruzados \(S_{xy}\), la cual a su vez se relaciona con la covarianza entre las variables \(x\) y \(y\).

Estos sencillos cálculos son realizados por R a través de la función “lm”, así:

# b1 = cov(x,y) / (sigma_x)²
b1 <- cov(x,y) / var(x)
print(b1)
## [1] 0.9449775
# bo = y_barra - b1*x_barra
b0 <- mean(y) - b1*mean(x)
print(b0)
## [1] 1436.679
residuales <- y - (b0 + b1*x)
n <- length(x)

sigma_2 <- sum(residuales^2)/(n-2)
print(sigma_2)
## [1] 288485.9
# Confirmamos usando la función lm, a continuación la documentación

# lm(formula, data, subset, weights, na.action,
#    method = "qr", model = TRUE, x = FALSE, y = FALSE, qr = TRUE,
#    singular.ok = TRUE, contrasts = NULL, offset, ...)

mod1 <- hass_dolar %>% lm(precio_aguacate_kg ~ precio_dolar,.)
print(mod1$coefficients)
##  (Intercept) precio_dolar 
## 1436.6790256    0.9449775
resumen <-  summary(mod1)
resumen$sigma
## [1] 537.1089
parametros <- data.frame(metodo = c("manual","lm"), 
                         beta_0 = c(b0,mod1$coefficients[1]),
                         beta_1 = c(b1,mod1$coefficients[2]),
                         sigma_2 = c(sigma_2, (resumen$sigma)^2))



parametros %>% 
  ftable()

metodo

beta_0

beta_1

sigma_2

manual

1,436.679

0.9449775

288,485.9

lm

1,436.679

0.9449775

288,485.9

#print(mod1$)
# Generando el modelo con R Base sería así:
# mod1 <- lm(promedio_aguacate_kg ~ promedio_dolar, hass_dolar)

Podemos notar entonces que ambos métodos son consistentes entre si, algo que no es de extrañar pues internamente R implementa estos métodos en sus librerías.

2. COMPROBACIÓN DEL SUPUESTO DE NORMALIDAD.

Hemos llegado al punto en que ya tenemos todos los elementos necesarios para definir formalmente el modelo estimado sobre nuestros datos. Podemos ahora concluir que el mejor modelo para describir la relación entre la variable precio del aguacate (\(y\)) y precio del dólar(\(x\)), pensado desde la perspectiva de minimización de cuadrados del error y estando restringidos a la familia de modelos de regresión lineal con componente de error normal es:

\(y = 1436.68 + 0.94x + e \text{ con e~N(0,288485.9)}\)

Pero, aunque nos hemos esforzado por encontrar las mejores estimaciones de los betas para minimizar la función de pérdida, esto no garantiza que nuestro modelo induzca errores normales, por tal motivo conviene recordar lo siguiente:

💭 Cuando planteamos el modelo, tuvimos que suponer algunas cosas, entre ellas, que la distribución de los errores sería normal. Por lo tanto, ahora que tenemos construído nuestro modelo, es necesario que validemos la suposición hecha.

🤔💭 ¿Cómo hacerlo ❔

Consideremos la siguiente tabla comparativa que informa sobre las ventajas y desventajas de las tres pruebas de normalidad (Kolmogorov-Smirnov, Lilliefors y Shapiro-Wilk) más usadas en estadística

tabla_normalidad <- data.frame(
  "Prueba" = c("Kolmogorov-Smirnov (KS)", "Prueba de Lilliefors", "Prueba de Shapiro-Wilk (SW)"),
  "Ventajas" = c(
    "- Puede utilizarse para cualquier distribución teórica.\n- Es una prueba no paramétrica versátil.\n - Adecuada para muestras grandes.",
    "- Es específica para evaluar la normalidad.\n- Puede ser más adecuada para muestras pequeñas.\n - Utiliza estadísticas de prueba que se ajustan a la distribución estándar.",
    "- Es especialmente potente para muestras pequeñas.\n- Diseñada específicamente para evaluar la normalidad.\n- Sensible a desviaciones de la normalidad en cualquier parte de la distribución."
  ),  "Desventajas" = c(
    "- Puede ser menos poderosa en muestras pequeñas.\n - No es específica para la distribución normal.",
    "- Restringida a la comprobación de normalidad. \n - Menos potente en muestras grandes.\n- No es tan versátil como KS para otras distribuciones.",
    "- Más compleja de entender y aplicar.\n- No proporciona estimaciones de parámetros.\n- Puede ser menos adecuada para muestras muy grandes.\n- No es tan versátil como KS para otras distribuciones."))

# Imprimir la tabla con formato

tabla_normalidad %>% 
  ftable() %>% 
  align(i = NULL, j = c(2:3),
        align = "left", part = "all")

Prueba

Ventajas

Desventajas

Kolmogorov-Smirnov (KS)

- Puede utilizarse para cualquier distribución teórica.
- Es una prueba no paramétrica versátil.
- Adecuada para muestras grandes.

- Puede ser menos poderosa en muestras pequeñas.
- No es específica para la distribución normal.

Prueba de Lilliefors

- Es específica para evaluar la normalidad.
- Puede ser más adecuada para muestras pequeñas.
- Utiliza estadísticas de prueba que se ajustan a la distribución estándar.

- Restringida a la comprobación de normalidad.
- Menos potente en muestras grandes.
- No es tan versátil como KS para otras distribuciones.

Prueba de Shapiro-Wilk (SW)

- Es especialmente potente para muestras pequeñas.
- Diseñada específicamente para evaluar la normalidad.
- Sensible a desviaciones de la normalidad en cualquier parte de la distribución.

- Más compleja de entender y aplicar.
- No proporciona estimaciones de parámetros.
- Puede ser menos adecuada para muestras muy grandes.
- No es tan versátil como KS para otras distribuciones.

Recordemos que la elección de la prueba depende de varios factores, incluyendo el tamaño de la muestra, el conocimiento previo sobre la distribución de los datos y los objetivos del análisis. Ninguna prueba es perfecta y es importante considerar el contexto y las características específicas de tus datos al seleccionar una prueba de normalidad.

🤔💭 Y si el texto anterior dice: “el conocimiento previo sobre la distribución de los datos” 🤔💭 ¿Será necesario representar la distribución en un gráfico? ¿Cómo funciona la prueba Lilliefors? A continuación describiremos los cálculos implicados:

# Paso 1: Ordenamos los datos de menor a mayor
n <- length(mod1$residuals)
df_residuals <- data.frame(residuals = mod1$residuals)
df_residuals <- df_residuals %>% 
                arrange(residuals)

# Paso 2: Estimamos los parámetros para la distribución normal teórica a partir de los datos

mean_residuals <- mean(df_residuals$residuals)
sd_residuals <- sd(df_residuals$residuals)

# Paso 3: Calculamos la función acumulada para cada valor de residual
cdf_teorica <- pnorm(df_residuals$residuals,
                      mean = mean_residuals,
                      sd = sd_residuals)

# Paso 4: Calculamos la máxima diferencia absoluta entre la densidad empírica y la teórica

# Método 1: Construir la función empírica usando fórmulas de R (ecdf) y luego tomar diferencias
#           absolutas entre valor empírico y teórico (enfoque KS)
cdf_empirica <- ecdf(df_residuals$residuals)
diferencias_absolutas <- abs(cdf_empirica(df_residuals$residuals) - cdf_teorica)
pos_estadistico_m1 <- which.max(diferencias_absolutas)
estadistico_m1 <- diferencias_absolutas[pos_estadistico_m1]

print(paste0("Este es el valor del estadístico: ", estadistico_m1," y esta es la posición en ",
             "la que ocurre la máxima diferencia ", pos_estadistico_m1))
## [1] "Este es el valor del estadístico: 0.0566890639519794 y esta es la posición en la que ocurre la máxima diferencia 60"
# Método 2: Construir manualmente los cálculos que permiten hallar la máxima diferencia
#           entre la empírica y la teórica usando la definición de la empírica como
#           función a trozos 📈

e_plus <- seq(1:n)/n
e_minus <- seq(0:(n-1))/n

D_plus <- e_plus - cdf_teorica
D_minus <- cdf_teorica - e_minus

pos_max_d_plus <- which.max(D_plus)
pos_max_d_minus <- which.max(D_minus)

estadistico_m2 <- max(D_plus[pos_max_d_plus], D_minus[pos_max_d_minus])
pos_estadistico_m2 <- function() {
  if(max(D_plus)>=max(D_minus)) {
    return(pos_max_d_plus)
    }else{
      return(pos_max_d_minus)
    }
  }

print(paste0("Este es el valor del estadístico: ", estadistico_m2," y esta es la posición en ",
             "la que ocurre la máxima diferencia ", pos_estadistico_m2()))
## [1] "Este es el valor del estadístico: 0.0566890639519794 y esta es la posición en la que ocurre la máxima diferencia 60"
# Comprobamos los resultados invocando la función de R que realiza estos cálculos automáticos

Lilliefors_test <- lillie.test(df_residuals$residuals)
print(Lilliefors_test)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  df_residuals$residuals
## D = 0.056689, p-value = 0.3637
estadistico_l <- Lilliefors_test$statistic
print(paste0("Este es el estadístico de Lilliefors usando paquetería de R: ", estadistico_l))
## [1] "Este es el estadístico de Lilliefors usando paquetería de R: 0.0566890639519794"

Hasta acá hemos verificado que los cálculos realizados manualmente coinciden con el código implementado por la librería nortest de R, y más específicamente en el código fuente de la función lillie.test. Ahora ilustremos algo de estos pasos usando una gráfica:

x_seg <- df_residuals$residuals[pos_estadistico_m1]
y_sup_seg <- cdf_teorica[pos_estadistico_m1]
y_inf_seg <- cdf_empirica(x_seg)

cat("x_seg: ", format(x_seg, digits=5), 
    " y_inf_seg: ", format(y_inf_seg, digits=5),
    " y_sup_seg: ", format(y_sup_seg, digits=5),
    sep='')
## x_seg: -147.97 y_inf_seg: 0.44776 y_sup_seg: 0.39107
datos_para_grafico <- data.frame(
  residuales = df_residuals$residuals
  )

ggplot(datos_para_grafico, aes(x = residuales)) +
#  geom_density(color = "blue", fill = "lightblue", alpha = 0.5) +
  stat_function(fun = pnorm, args = list(mean = mean_residuals,
                                         sd = sd_residuals),
                color = "red", size = 1) +
  stat_ecdf(color = "green", size = 1) +
  geom_segment(aes(x = x_seg, y = y_inf_seg, xend = x_seg,
                   yend = y_sup_seg),
               color = "purple", size = 1, linetype = "solid") +
  geom_point(aes(x = x_seg, y = y_inf_seg),
             color = "purple",size = 3) +
  geom_point(aes(x = x_seg, y = y_sup_seg),
             color = "purple", size = 3) +
  annotate("text", x = x_seg , y = (y_inf_seg + y_sup_seg)/2, 
           label = paste0("Máxima Separación=",
                          format(estadistico_m1,
                                 digits=5)," 🤔💭", sep=''),
           color = "purple", size = 3, hjust = -0.05, vjust = 0) +
  coord_cartesian(xlim = c(min(datos_para_grafico$residuales),
                           max(datos_para_grafico$residuales))) +
  xlab("Residuales del modelo") +
  ylab("Probabilidad acumulada") +
  labs(title = "Esquema de funcionamiento de la prueba de Lilliefors") +
  theme_minimal()

Nos queda la siguiente reflexión sobre la gráfica anterior, ¿Por qué elegir la máxima diferencia absoluta entre las dos curvas en lugar de un promedio de las diferencias a lo largo de todo el eje? ¿Por qué no usar una media ponderada de los residuales al cuadrado?

Comparemos los resultados contra la prueba de Shapiro WIlk y Kolmogorov

# Probando normalidad del error. 

shapiro.test(mod1$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  mod1$residuals
## W = 0.98361, p-value = 0.1078
ks.test(mod1$residuals, "pnorm")
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  mod1$residuals
## D = 0.50746, p-value < 0.00000000000000022
## alternative hypothesis: two-sided

Hemos notado diferencias con la prueba de Kolmogorov, ¿Cuál es la razón?, podrás encontrar en el siguiente artículo en el que Lilliefors introduce por primera vez su prueba:

Lilliefors, H. (June 1967), “On the Kolmogorov–Smirnov test for normality with mean and variance unknown”, Journal of the American Statistical Association, Vol. 62. pp. 399–402.

Una aclaración respecto a que esta prueba es menos potente ❗ en comparación con la prueba KS estándar, esto significa que es menos sensible para detectar desviaciones de la normalidad en muestras grandes, especialmente si los datos no siguen una distribución estrictamente normal. La prueba KS estándar es más robusta en este sentido, pero por lo mismo más conservador. En el contexto de modelos de regresión, no requerimos un cumplimiento estricto de la normalidad sino una semejanza razonable para habilitar el uso del modelo. En ese sentido es preferible usar Lilliefors.😀😎

🤔💭 Otras posibles preguntas son:

  • ¿Cuáles son los riesgos que implica asumir a priori un comportamiento específico para el error?
  • ¿Si el modelo no cumple con lo supuesto para el error, qué alternativas tengo?
  • Al plantear un modelo alternativo que captura mejor el comportamiento del error usando una distribución distinta a la normal ¿Puedo seguir usando los métodos de inferencia tradicionales o qué modificaciones debo introducir?
  • ¿Por qué es importante para un magister en estadística conocer sobre métodos de simulación Montecarlo?

3. EL SUPUESTO DE INDEPENDENCIA DEL ERROR

Un supuesto que hemos perdido de vista por un momento en los apartados anteriores es el supuesto de independencia del error. Su consideración es importante porque representa un punto de quiebre en el uso de modelos de regresión vs modelos alternativos como los de series de tiempo. Pero ¿Qué significa exactamente la independencia?. Consideremos los siguientes escenarios:

  • Un modelo que presenta sistemáticamente más dificultades para predecir de forma exacta valores altos de la variable respuesta pero es bueno prediciendo valores bajos.

  • Un modelo que se construyó sobre datos de un experimento de durezas de metales, donde el mecanismo que sensa la dureza perdió calibración o exactitud después de tomar la medida de las 100 primeras muestras.

  • Un modelo que intento predecir un precio de activo financiero en función de un activo de referencia, pero conforme pasa el tiempo el modelo se va desajustando y haciendo menos preciso.

  • Un modelo que desconce la presencia de una variable oculta que controla la variación acoplada de mi variable respuesta vs mi variable predictora y las hace parecer más correlacionadas de lo que realmente son (regresiones espurias)

Para cada uno de los casos anteriores el investigador tiene que decidir el tipo de recurso gráfico o estadístco que le permita inspeccionar si una situación de esa naturaleza está ocurriendo. Las alternativas son:

Gráficos de residuales vs. valores ajustados: Graficar los residuales (diferencia entre los valores observados y los valores predichos) contra los valores ajustados (las predicciones del modelo) es una forma común de detectar patrones de no independencia.

🤔💭 ¿Qué comportamiento es deseable para este gráfico?

Si no hay patrones evidentes en el gráfico y los residuales se distribuyen aleatoriamente alrededor de cero, es un indicio de independencia.

Gráfico de residuales con índice temporal: Si tus datos están ordenados en el tiempo (series temporales), un gráfico de residuales en el tiempo puede revelar patrones de autocorrelación. Deben parecer ruido blanco, es decir, sin patrones visibles. También es conveniente usar gráficos de autocorrelación de residuales (ACF) y gráficos de autocorrelación parcial de residuales (PACF), ya que estos informan sobre una posible correlación entre el error actualmente observado y errores observados en otros momentos del pasado.

Prueba de Durbin-Watson: Esta prueba estadística evalúa si existe autocorrelación de primer orden en los residuales, mejora el análisis gráfico por ser una prueba formal. Un valor de Durbin-Watson cercano a 2 indica independencia, mientras que valores significativamente diferentes de 2 pueden indicar autocorrelación. De nuevo es necesario elegir una variable de ordenamiento para los errores encaminada a detectar el tipo de patrón de interés.

🤔💭 ¿Cómo hacer esto?

  • Ordene los errores según magnitud del valor predicho (ajustado)
  • Ordene los errores según el tiempo.
  • Ordene los errores según el orden de experimentación.
  • No ordene por la posición del registro en el data frame a menos qué entienda qué significa ese orden y qué aporta al análisis.
  • Ordene por una variable de interés en la que quiera inspeccionar si los residuales se ven afectados por la magnitud de tal variable.

👦👴👧👨👨‍👵

Ejemplo: Puede considerar revisar si el consumo calórico puede explicar residuales de un modelo simplificado en el que la estatura pretenda explicar el peso. 🔎🕵️ Aunque este ejemplo parece obvio ilustra una idea muy importante, porque nos da un criterio para ingresar una nueva variable a un modelo de regresión: Ingrésela cuando dicha variable se relaciona bien con los residuales generados por la primera versión de mi modelo!!!. Es decir, la variable debe ingresar al modelo cuando pueda explicar la variabilidad restante no explicada por mi modelo anterior.

Prueba de Ljung-Box: Esta prueba se utiliza para evaluar la autocorrelación en los residuales en varios retardos. Si los residuales son independientes, las autocorrelaciones en los retardos deberían ser cercanas a cero. Por varios retardos entendemos inspeccionar relación no sólo con el valor del residual inmediatamente anterior, sino contra cualquier otro ubicado en un momento más anterior. 🔎 Es una mejora a la prueba de Durbin-Watson en la medida en que no limita la autocorrelación al retardo anterior.

Estas son algunas de las formas comunes de comprobar la independencia de los errores en un modelo de regresión lineal. Si encuentras evidencia de autocorrelación o patrones en los residuales, puede ser necesario considerar modelos más avanzados, como 📉 🚪modelos autorregresivos, o de media móvil, o en forma más general modelos de series temporales, para capturar adecuadamente la estructura de dependencia en los datos.

En conclusión:

En regresión lineal asumimos la no correlación de los errores. En nuestro caso el precio del aguacate y el precio del dólar son variables evolucionando en el tiempo, así que quisiéramos que tras el transcurrir del tiempo nuestro modelo no se deteriore. Una forma en que esto puede ocurrir es que alguna autocorrelación temporal provoque residuales que van variando en sus propiedades conforme el tiempo pasa. Por tal motivo necesitamos técnicas como las anteriores para probar independencia.

Por ejemplo puede ocurrir que el error se vaya volviendo más grande en promedio con el paso del tiempo o vaya creciendo en dispersión, o cualquier otra anomalía que impida suponer que tenemos un modelo estable para todo momento del tiempo. Los gráficos de dispersión típicos ilustran la asociación entre la magnitud de la variable predictora y la variable respuesta, pero para probar independencia es preferible que los errores se dibujen indexados por la fecha en la que tal error ocurrió, y así poder apreciar posibles patrones temporales. Haremos esas inspecciones a continuación 🔎 :

# Primero inspeccionamos las variables individuales para observar su
# comportamiento temporal en un mismo gráfico

# library(lubridate)
x <- hass_dolar$precio_dolar
y <- hass_dolar$precio_aguacate_kg

hass_dolar$fecha <-as.Date(paste0(hass_dolar$anio,"-",hass_dolar$mes,"-01"))
hass_dolar$predicciones <- mod1$fitted.values
hass_dolar$residuales <-mod1$residuals
  
# Gráfico de evolución de precios y predicciones

ggplot(hass_dolar, aes(x = fecha)) +
  geom_line(aes(y = precio_dolar, color = "Precio del dólar"), size = 1) +
  geom_line(aes(y = precio_aguacate_kg, color = "Precio del aguacate"), size = 1) +
  geom_line(aes(y = predicciones, color = "Predicciones precio aguacate"), size = 1) +
  scale_color_manual(values = c("Precio del dólar" = "blue",
                                "Precio del aguacate" = "green",
                                "Predicciones precio aguacate" = "red")) +
  labs(title = "Evolución temporal del precios y predicciones",
       x = "Fecha",
       y = "Valor") +
  theme_minimal() +
  scale_x_date(date_breaks = "1 year", date_labels = "%b %Y") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

De entrada, podemos observar en el gráfico anterior que hay zonas donde el precio del aguacate fue más volátil que el precio del dólar, razón por la cual el precio del dólar no pudo seguir eficientemente el patrón del precio del aguacate, veremos esto cómo afectó los residuales:

# Gráfico de residuales vs valores ajustados
hass_dolar %>% 
  ggplot(aes(x= predicciones, y=residuales)) +
  geom_point()+ 
  theme_light()+
  ylab("Residuales") +
  xlab("Valor predicho")

# Gráfico de residuales con índice temporal
hass_dolar %>% 
  ggplot(aes(x= fecha, y=residuales)) +
  geom_point()+ 
  geom_line(aes(y = residuales, color = "residuales"), size = 1)+
  theme_light()+
  xlab("Fecha")+
  ylab("Residuales")

# Prueba de Durbin-Watson tomando el mismo orden en el que venían los
# datos en el dataframe puesto que este orden se corresponde con el 
# orden temporal.

resultado_dw <- dwtest(mod1)
print(resultado_dw)
## 
##  Durbin-Watson test
## 
## data:  mod1
## DW = 0.47441, p-value < 0.00000000000000022
## alternative hypothesis: true autocorrelation is greater than 0
# Prueba de Ljung-Box
# Como esperamos cierta estacionalidad intra año en los datos podría tener sentido esperarla ttambién para los residuales. Por lo tanto eligiremos lag = 1,2,3,...,12

p_values <- vector()
for (i in 1:12) {
  resultado_ljung_box <- Box.test(mod1$residuals, lag = i,
                                  type = "Ljung-Box")
  p_values[i] <-resultado_ljung_box$p.value
}

resultados_lb <- data.frame(
  lag = seq(1,12),
  p_value = p_values
)

print(resultado_ljung_box)
## 
##  Box-Ljung test
## 
## data:  mod1$residuals
## X-squared = 124.84, df = 12, p-value < 0.00000000000000022
print(resultados_lb)
##    lag p_value
## 1    1       0
## 2    2       0
## 3    3       0
## 4    4       0
## 5    5       0
## 6    6       0
## 7    7       0
## 8    8       0
## 9    9       0
## 10  10       0
## 11  11       0
## 12  12       0

Conclusiones:

  • Los errores no presentan un patrón claro al momento de compararse contra los valores ajustados, es decir que el modelo tiene una calidad similar prediciendo valores pequeños, intermedios y grandes en el precio del aguacate.

  • El gráfico de residuales indexados por el tiempo parece tener ciertas oscilaciones predecibles pero es difícil de establecer por simple inspección visual.

  • La prueba de Durbin Watson nos ayuda a revisar si hay autocorrelación del error de hoy con el error de hace un mes (recuerde que la serie es mensual). De hecho la prueba arroja un estadístico de prueba \(DW < 2\) indicando con ello autocorrelación positiva de los errores.

  • En consistencia con lo anterior la prueba de Ljung-Box arroja valores p sistemáticamente menores que 0.05 para todos los rezagos considerados: 1 mes, 2 meses, 3 meses, …, 1 año. Y decimos que en consistencia con la prueba de Durbin Watson porque la prueba de Ljung box será indicativa de autocorrelación en máximo n rezagos si ya lo es para máximo n=1 rezagos como lo indica la prueba DW.

  • La no independencia del error invalida completamente mi modelo, porque en estos escenarios la varianza del error del modelo se encuentra subestimada, pues ante la presencia de autocorrelación de los residuales no podemos afirmar que \(\sigma_e² = var(error)\) ya que esto asume implícitamente independencia. En procesos autocorrelacionados la varianza es realmente mayor que esto. Usar un error subestimado puede llevarme a conclusiones falsamente significativas, como por ejemplo puede declarar la variable predictora como significativa, al estimar inadecuadamente el error para \(\beta_1\) ya que \(\widetilde{V}(\widetilde{\beta}_1)\) depende de \(\widetilde{\sigma}²_e\)

Dejaremos el estudio de las cuestiones relativas al tratamiento de errores autocorrelacionados para el capítulo de series de tiempo.

4. EL SUPUESTO DE HOMOCEDASTICIDAD

De manera informal este supuesto busca verificar que la varianza de los errores es estable, bien sea respecto al tiempo, respecto a una variable predictora, o respecto a la magnitud de los valores predichos. Las inspecciones de estas cuestiones usualmente se basan en el estudio de los patrones de los gráficos de residuales, en específico estudiamos:

  • Si un gráfico de residuales estandarizados (divididos por su desviación estándar) vs valor predicho arroja dispersiones diferentes conforme el valor predicho crece. En un gráfico como éste si observamos que la dispersión de los puntos cambia a lo largo de los valores ajustados (por ejemplo, los puntos se abren o se estrechan a medida que aumentan los valores ajustados), eso podría indicar la presencia de heterocedasticidad.

  • Si el modelo posee una variable categórica predictiva (no es nuestro caso), respecto a la cual podamos validar si la dispersión del error es mayor en los grupos definidos por esa variable predictiva. Ampliaremos estas cuestiones en el apartado de modelo de regresión lineal múltiple.

  • Si un gráfico de residuos vs variables independientes revela que la dispersión de los residuos cambia con respecto a cada variable independiente. Por ejemplo, podemos patrones de abanico o embudo que puedan indicar problemas de heterocedasticidad.

  • Si alguna pruebas estadísticas formales como la prueba de Breusch-Pagan o la prueba de White proporcionan evidencia de heterocedasticidad al comprobarse que los cuadrados de los residuales se relación sistemáticamente con el conjunto de variables predictoras, esto es, si es posible construir un modelo no ingenuo que haga depender el cuadrado del error de los valores de las variables predictoras. Si tal cosa existe significa que los errores son heterocedasticos en la medida en que su varianza depende o cambia según los valores de las variables predictoras. Apliquemos estos análisis a continuación:

# Residuales estandarizados vs valores predichos
hass_dolar %>% 
  ggplot(aes(x= predicciones, y=residuales/sd(residuales))) +
  geom_point()+ 
  theme_light()+
  ylab("Residuales estandarizados") +
  xlab("Valor predicho")

# Residuales estandarizados vs precio del dólar
hass_dolar %>% 
  ggplot(aes(x= precio_dolar, y=residuales/sd(residuales))) +
  geom_point()+ 
  theme_light()+
  ylab("Residuales estandarizados") +
  ylab("Precio del dólar")

Como se puede apreciar, el gráfico de errores estandarizados vs valores predichos será muy parecido al gráfico de errores vs variable predictora en un modelo de regresión lineal simple, puesto que en este contexto las predicciones dependen de la variable de pronóstico en forma tal que la primera no es más que un rescalado y desplazamiento de la segunda. Tendría más sentido en un modelo de regresión lineal múltiple generar ambos gráficos, acá uno sólo de los anteriores bastaría.

Por otro lado, ambos gráficos revelan que no existe un patrón específico que haga sospechar de la presencia de heterocedasticidad, hagamos un inspección más cuidadosa usando pruebas formales:

# Pruebas estadísticas formales
resultado_breusch_pagan <- bptest(mod1)
print(resultado_breusch_pagan)
## 
##  studentized Breusch-Pagan test
## 
## data:  mod1
## BP = 1.741, df = 1, p-value = 0.187

La prueba de Breusch-Pagan propuesta en (1979) consiste en ajustar un modelo de regresión lineal con variable respuesta dada por residuales del modelo original al cuadrado \(e_i²\) y como covariables las variables del modelo original. Por ejemplo, si se tienen \(k=2\) covariables para explicar a \(Y\) entonces el modelo de regresión para estudiar la homocedasticidad es:

\[e_i² = \delta_0 + \delta_1*x_1 + \delta_2*x_2 + u\]

Si se concluye que \(\delta_1 = \delta_2 = 0\), significa que los residuales no son función de las covariables del modelo. El estadístico en esta prueba está dado por \(n*R²\) y bajo la hipótesis nula verdadera, el estadístico tiene distribución \(\chi²_k\)

Lo que el resultado nos está indicando es que \(\chi²_1 = 1.741\), y \(df = 1\), con \(\text{p value} = 0.187\), y debido a que p>0.05, no hay evidencia estadística suficiente para el rechazo de la hipótesis de homocedasticidad. Todo esto en concordancia con lo revelado por el análisis gráfico.

Sin embargo, el test de Breusch-Pagan sólo detecta formas lineales de heterocedasticidad, esto significa que sólo relaciona los cuadrados del error con términos que dependen linealmente de las covariables. Para resolver esta limitación, el test de White propuesto por White (1980), permite contrastar no linealidades utilizando los cuadrados y los productos cruzados de todos los regresores. Sea que se use Breusch-Pagan o White, note que sólo se están explorando dos formas de heterocedasticidad en un universo infinitamente amplio de formas en que ésta se puede presentar, esto es, sólo se explora heterocedasticidad con forma parabólica y con forma lineal. Se podrían explorar otras pero conforme la complejidad del modelo crece se pierde su utilidad en la práctica.

En el caso del test de White, si \(k=2\) el modelo que relaciona \(e_i²\) con las covariables queda así:

\[e_i² = \delta_0 + \delta_1x_1 + \delta_2x_2 + \delta_3x_1x_2 + \delta_4 x_1² + \delta_5x_2² + u\] Podemos pensar el modelo de White como una extensión del modelo de Breusch-Pagan, por eso no es de extrañar que en la siguiente celda usemos la misma función en R, pero agregando un término al cuadrado para la única covariable implicada en nuestro ejemplo. No habrán términos cruzados porque sólo tenemos una variable predictora.

resultado_white <- bptest(mod1,
                          varformula = ~ (precio_dolar +
                                            I(precio_dolar^2)),
                          data=hass_dolar)
print(resultado_white)
## 
##  studentized Breusch-Pagan test
## 
## data:  mod1
## BP = 7.8925, df = 2, p-value = 0.01933

Para ilustrar como trabajan estas pruebas por debajo, considere la siguiente implementación manual:

# Implementación manual del test de White

#1. Construir el dataframe que relaciona el cuadrado del error con
#los términos lineales de la predictora y el término cuadrático
#(por ser modelo lineal simple no hay productos cruzados)

df_squared_error <- data.frame(
  squared_error = mod1$residuals ^ 2,
  precio_dolar = hass_dolar$precio_dolar,
  precio_dolar_2 = hass_dolar$precio_dolar^2
)

#2. Calcule un modelo lineal ajustado sobre el cuadrado del error regresado por el termino lineal y cuadrático del precio del dolar

modelo_white <- lm(data = df_squared_error,
                   squared_error ~ precio_dolar + precio_dolar_2)

#3. Calcule el estadístico de prueba como n*R^2 (con n=datos del dataframe de errores)
R2 <- summary(modelo_white)$`r.squared`
w_estadistico = nrow(df_squared_error)*R2
#4. Como el estadístico se distribuye ji-cuadrado con 2 grados de libertad (2 variables regresoras) se tiene que el p valor es:
p_valor = 1-pchisq(as.numeric(w_estadistico),2)
paste("BP = ",round(w_estadistico,4)," y p_valor=",round(p_valor,5))
## [1] "BP =  7.8925  y p_valor= 0.01933"

En la salida impresa puede notar que los resultados son consistentes con lo generado por la prueba de heterocedasticidad arrojada con la función bptest.

Además, cuando se estudia heterocedasticidad en segundo grado se aprecia que el modelo tiene varianzas del error no constantes porque estas dependen del cuadrado del precio del dólar. Veamos un poco esto:

estimated_squared_error <- modelo_white$coefficients[1] + modelo_white$coefficients[2]*df_squared_error$precio_dolar +modelo_white$coefficients[3]*df_squared_error$precio_dolar_2

p<- ggplot(data = df_squared_error, aes(x=precio_dolar,
                                        y=squared_error))+
  geom_point()+
  geom_line(aes(y=estimated_squared_error, color="red"))+
  labs(title="Modelo white sobre cuadrados del error",
         subtitle = paste("R²=",round(R2,4),
                          "W-Estadístico=",round(w_estadistico,4),
                          "p_valor=",round(p_valor,5)))+
  scale_color_manual(values=c("red"="blue"),
                     labels=c("Error cuadrado estimado"),
                     name="Modelo de White")

print(p)

En este caso \(W = nR² = 7.8925\) es el estadístico de la prueba con \(k=2\) grados de libertad asociados a el número de covariables. Ahora bien \(\text{p valor} = 0.01933\), obtenido a través de calcular \(p(\chi^2_2>=7.8925)\), indicando esta vez que hay evidencia estadística para asumir heterocedasticidad de orden cuadrático. Podrá notar que un patrón parabólico no muy fuerte se observa cuando dibujamos el modelo sobre la nube de puntos del error cuadrado vs el precio del dólar.

Vamos a intentar resolver el problema de heterocedasticidad con una transformación de la variable respuesta que tienda a aplanar la curva de errores cuadrados responsable de la falla del supuesto. Dos posibles transformaciones son:

    1. Extraer raíz cuadrada de la respuesta.
    1. Extraer logaritmo natural de la respuesta.

La segunda transformación es mucho más severa. Probaremos con esta.

# Datos para nuevo modelo con y modificada a través de logaritmo natural
data_mod <- data.frame(
  y_mod = log(hass_dolar$precio_aguacate_kg),
  precio_dolar = hass_dolar$precio_dolar,
  precio_dolar_2 = hass_dolar$precio_dolar^2
)

# Nuevo modelo lineal
mod2<-lm(y_mod~precio_dolar, data=data_mod)
residuales_mod<-mod2$residuals

#Guardamos los cuadrados del nuevo error
data_mod$squared_error <- residuales_mod^2

#Aplicamos prueba de white sobre el modelo con la y modificada
white_test_mod<-bptest(mod2,
                       varformula = ~ (precio_dolar + I(precio_dolar^2)),
                       data=data_mod)
print(white_test_mod)
## 
##  studentized Breusch-Pagan test
## 
## data:  mod2
## BP = 5.2735, df = 2, p-value = 0.07159
w_estadistico_mod <- white_test_mod$statistic
p_value_mod <- white_test_mod$p.value

#Obtenemos los coeficientes para dibujar curva de errores cuadrados
#estimados

manual_white_test_mod <-lm(data = data_mod,
                           squared_error ~ precio_dolar + precio_dolar_2)

coefs_mod <- manual_white_test_mod$coefficients
R2_mod <- summary(manual_white_test_mod)$`r.squared`

#Obtenemos los nuevos errores cuadrados estimados
data_mod$estimated_squared_error <- (coefs_mod[1] +
                                    coefs_mod[2]*data_mod$precio_dolar +
                                    coefs_mod[3]*data_mod$precio_dolar^2)

#Comprobamos que no se nos deteriore la normalidad
Lilliefors_test_mod <- lillie.test(mod2$residuals)
print(Lilliefors_test_mod)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  mod2$residuals
## D = 0.050708, p-value = 0.544
#Comprobamos si milagrosamente resolvimos la no independencia del modelo
#original
resultado_ljung_box_mod <- Box.test(mod2$residuals, lag = i,
                                    type = "Ljung-Box")
print(resultado_ljung_box_mod)
## 
##  Box-Ljung test
## 
## data:  mod2$residuals
## X-squared = 133.25, df = 12, p-value < 0.00000000000000022
#Graficamos de nuevo el modelo de white sobre los cuadrados del error
#para el modelo con y modificada
p<- ggplot(data = data_mod, aes(x=precio_dolar,
                                  y=squared_error))+
  geom_point()+
  geom_line(aes(y=estimated_squared_error, color="red"))+
  labs(title="Modelo white sobre cuadrados del error",
         subtitle = paste("R²=",round(R2_mod,4),
                          "W-Estadístico=",round(w_estadistico_mod,4),
                          "p_valor=",round(p_value_mod,5)))+
  scale_color_manual(values=c("red"="blue"),
                     labels=c("Error cuadrado estimado"),
                     name="Modelo de White")

print(p)

Las conclusiones de todo esto son:

  • Nuestro modelo es inservible, porque la independencia sigue sin cumplirse, supuesto fundamental para poder usar el modelo con la confianza de no estar subestimando la varianza del error.
  • La prueba de homocedasticidad pasa por un estrecho margen pese a las severas transformaciones.
  • No tenemos problemas con normalidad.

Hasta acá las discusiones del modelo lineal simple con este dataset de ejemplo. Para nuestro capítulo siguiente relacionado con pruebas de hipótesis para los betas del modelo, tendremos que usar un modelo que cumpla con todos los supuestos. Para este efecto aprovecharemos para simular datos bajo un modelo teórico perfecto. Esta herramienta de simulación permitirá tomar control del cumplimiento de los supuestos para investigar tranquilamente el comportamiento estadístico de los betas estimados de un modelo cuando todos las condiciones para su uso están garantizadas.

5. PRIMER TALLER

  1. Calcular un modelo de regresión para los cambios porcentuales del precio del aguacate de un mes a otro vs los cambios porcentuales del precio del dólar de un mes a otro. ¿En qué difiere este enfoque del originalmente seguido? ¿Cuál cree que sea la ventaja de generar un modelo de esta naturaleza?

  2. Evalue la normalidad del error para este modelo usando la prueba de Lilliefors

  3. Evalúe en un gráfico la consistencia de la prueba de Lilliefors (Pista: Use el código ya suministrado para calcular la máxima distancia entre la curva normal teórica y la curva empírica)

  4. Pruebe la homocedasticidad del nuevo modelo usando tanto gráficos como pruebas formales, y dibuje la curva de errores cuadrados estimados basada en el modelo de Breusch-Pagan, es decir, sin incluir el término cuadrático. ¿Es esta prueba suficientemente concluyente?

  5. Pruebe la independencia de los errores del modelo, comparéla contra los resultados obtenidos con el modelo que usaba los precios en valores absolutos, es decir, el modelo discutido en este documento.

  6. Extraiga conclusiones generales de la adecuación del nuevo modelo a todos los supuestos, ¿Cuál de los dos modelos sería más útil para hacer predicciones, el modelo elaborado por los profesores o el que ustedes elaboraron? Felicítese mucho si el suyo es mejor. :)

6. PRUEBAS DE HIPÓTESIS SOBRE LOS BETAS DEL MODELO

Para poder entender muy bien lo que significa probar hipótesis sobre los betas de un modelo de regresión, lo primero que hay que tener en claro es que un modelo teórico se puede pensar como una máquina generadora de datos 🧮. Como en toda máquina habrá una ajuste inicial que requiere ser definido antes de producir algo con la máquina, pero el problema radica en que usualmente no estamos presentes en el momento en que tal máquina se ajusta, sino que debemos inferir 🕵️‍♂️el ajuste de acuerdo con los datos que la máquina ya nos está ya arrojando.

Es decir, debemos reconocer la incertidumbre que tenemos frente al conocimiento del ajuste inicial e intentar minimizarla con base en los datos ya generados. Otra cuestión es que no podemos quedarnos infinitamente observando los datos que la máquina genera, sino que usualmente tenemos una muestra limitada de los mismos. Las preguntas centrales que el estadístico debe responder son: 🤔

  • ¿Cuáles son estas configuraciones iniciales en el caso de un modelo de regresión?
  • ¿Cómo podemos aprovechar al máximo la muestra de datos observados para intentar inferir la configuración inicial de la máquina (de ahora en adelante llamada más formalmente proceso generador de datos)?
  • ¿Qué propiedades deseables debería tener la técnica que usemos para aprovechar la información limitada de la muestra con la que contamos para inferir la configuración inicial(de ahora en adelante más formalmente llamada parámetros)?

Retomemos entonces nuestros modelo, para aterrizar mejor estas cuestiones:

\[Y = \beta_0 + \beta_1X + e\] Nuestro proceso generador de datos consiste así en una función lineal que hace depender a \(Y\) del valor de \(X\) y de un componente aleatorio \(e\) que indica la parte del valor de Y que no logra ser capturada por la relación lineal que la X tiene con la Y. Nosotros no conocemos \(\beta_0\) y \(\beta_1\), así como tampoco \(\sigma\) (recordemos que \(e\text{~}N(0,\sigma²)\)). Así que debemos usar la información de la muestra de datos observados \((x_1, y_1), (x_2, y_2), ..., (x_n, y_n)\) para intentar deducir los valores de estos tres parámetros. Sin embargo tales valores son sólo un estimado de los verdaderos parámetros del proceso generador de datos.

En capítulos anteriores dejamos claro que un método posible para estimar los parámetros \(\beta_1\) y \(\beta_2\) es usar minimización del MSE(Mean Squared Error), es un esfuerzo por generar la mejor recta ajustada a los datos, pero cabe preguntarse:

  • ¿Generar una recta que se ajuste lo mejor posible a la muestra garantiza obtener estimaciones que realmente se parezcan a los parámetros de mi proceso generador de datos?

Lo que queremos indicar con esta pregunta es que: un buen ajuste no tiene que ser sinónimo de buena estimación de parámetros, y si lo es deberíamos poderlo demostrar, porque el ajuste habla de la calidad del modelo para explicar los datos observados, y no necesariamente para explicar todo el universo de los datos que mi proceso generador de datos podría generar. Si nos ajustamos bien a la muestra, esto no implica a-priori ajustarnos bien a la población. No obstante algo nos hace sospechar que no es tan mala idea usar una estrategia como esta porque después de todo no conocemos más datos que los que la muestra pueda revelarnos.

6.1 Estudiando el comportamiento estadístico de los estimadores

¿Cómo garantizamos que dentro de todas las posibles operaciones matemáticas que podemos ejecutar con los valores de la muestra, las operaciones:

\[\widetilde{\beta}_1 = \frac{S_{xy}}{S_{xx}}\] \[\widetilde{\beta}_0 = \bar{y} - \widetilde{\beta}_1\bar{x}\] Sean las operaciones que conducen a valores más próximos a \(\beta_1\) y \(\beta_0\) respectivamente, aun si la muestra de datos cambia?.

No podemos demostrar esto con el rigor necesario porque escapa un poco del alcance del presente curso, pero podemos al menos ilustrar lo que queremos indicar por cercanía a los betas del modelo usando el siguiente enfoque analítico.

Detengamonos un instante a analizar la estructura de \(\widetilde{\beta}_1\), notemos que este estimador puede pensarse como una combinación lineal de las \(y_i\), veamos por que:

\[\widetilde{\beta}_1 = \frac{S_{xy}}{S_{xx}}=\frac{\sum_{i=1}^{n}(x_i-\bar{x})*(y_i-\bar{y})}{S_{xx}}=\sum_{i=1}^{n}\frac{y_i(x_i-\bar{x})-\bar{y}(x_i-\bar{x})}{S_{xx}}\\=\sum_{i=1}^{n}\frac{y_i(x_i-\bar{x})}{S_{xx}}-\sum_{i=1}^{n}\frac{\bar{y}(x_i-\bar{x})}{S_{xx}}\\=\sum_{i=1}^{n}\frac{y_i(x_i-\bar{x})}{S_{xx}}-\frac{\bar{y}}{S_{xx}}\sum_{i=1}^{n}(x_i-\bar{x})\\=\sum_{i=1}^{n}\frac{y_i(x_i-\bar{x})}{S_{xx}}-\frac{\bar{y}}{S_{xx}}\left(\sum_{i=1}^{n}x_i-\sum_{i=1}^{n}\bar{x}\right)\\=\sum_{i=1}^{n}\frac{y_i(x_i-\bar{x})}{S_{xx}}-\frac{\bar{y}}{S_{xx}}\left(n\bar{x}-n\bar{x}\right)=\sum_{i=1}^{n}\frac{y_i(x_i-\bar{x})}{S_{xx}}\] Ahora sea \(c_i= \frac{x_i-\bar{x}}{S_{xx}}\) entonces concluimos que:

\[\widetilde{\beta}_1=\sum_{i=1}^{n}c_i*y_i\] Es decir \(\widetilde{\beta}_1\) es en efecto una combinación lineal de los \(y_i\) con coeficientes \(c_i\), por lo tanto, retomando a nuestro proceso generador de datos, y tomando un \(y_i\) particular generado por dicho proceso, tenemos que \(y_i = \beta_0 + \beta_1x_i + e_i\) donde los \(e_i\) se asumen normales e independientes, lo cual implica que:

\[E(y_i|x_i) = E((\beta_0 + \beta_1x_i + e_i)|x_i)=\beta_0+\beta_1x_i + E(e_i) = \beta_0 + \beta_1x_i\]

Ya que \(E(e_i)=0\) desde los supuestos. Y por otro lado que:

\[Var(y_i|x_i) = Var((\beta_0 + \beta_1x_i + e_i)|x_i) = 0+0+Var(e_i)= \sigma^2\]

Ya que \(Var(e_i) = \sigma^2\) también desde los supuestos.

Así que como conclusión tenemos que, cuando el \(x_i\) está ya dado:

\[y_i\text{~}N(\beta_0 + \beta_1x_i,\sigma^2)\]

Todo esto por propiedad de linealidad de la normalidad, es decir, debido a que una transformación afín de la variable normal \(e_i\) dada por un desplazmiento \(\beta_0\) y un factor de escala \(\beta_1\) es también una variable aleatoria normal.

Lo anterior, tácitamente nos está diciendo que nuestro estimador \(\widetilde{\beta}_1\) es entonces una combinación lineal de variables aleatorias normales, razón por la cual, partiendo de la propiedad reproductiva de la normal podremos afirmar que, cuando todos los \(x_i\) están dados \(\widetilde{\beta}_1\) es también una variable aleatoria normal.

Investiguemos su media:

\[E(\widetilde{\beta}_1)=E\left(\sum_{i=1}^{n}c_iy_i\right)=\sum_{i=1}^{n}E(c_iy_i)=\sum_{i=1}^{n}c_iE(y_i)=\sum_{i=1}^{n}c_i(\beta_0 + \beta_1x_i)\\=\beta_0\sum_{i=1}^{n}c_i + \beta_1\sum_{i=1}^{n}c_ix_i\] Pero es fácil ver que \(\sum_{i=1}^{n}c_i=0\) y que \(\sum_{i=1}^{n}c_ix_i=1\) (por favor demuéstrelo), y entonces:

\[E(\widetilde{\beta}_1)=\beta_1\]

Lo que indicaría que \(\widetilde{\beta}_1\) es un estimador insesgado para \(\beta_1\)!!!

Y ahora investiguemos su varianza:

\[Var(\widetilde{\beta}_1)=Var\left(\sum_{i=1}^{n}c_iy_i\right)=\sum_{i=1}^{n}{c_i}²Var(y_i)=\sum_{i=1}^{n}{c_i}²{\sigma}²={\sigma}²\sum_{i=1}^{n}{c_i}²\] Ya que de nuevo por asunciones sobre el proceso generador de datos las \(y_i\) son no correlacionadas y habíamos demostrado previamente que \(Var(y_i) = {\sigma}²\), sin importar el \(y_i\) considerado, es decir, estos son homocedásticos como consecuencia de la definición del modelo. Por lo tanto:

\[Var(\widetilde{\beta}_1)={\sigma}²\sum_{i=1}^{n}{\left(\frac{x_i-\bar{x}}{S_{xx}}\right)}^2=\frac{\sigma²}{{S²}_{xx}}\sum_{i=1}^{n}(x_i-\bar{x})²=\frac{\sigma²}{{S²}_{xx}}*S_{xx}=\frac{\sigma²}{{S}_{xx}}\] Finalmente, ya que hemos investigado ampliamente el comportamiento del estimador \(\widetilde{\beta}_1\) usemos este conocimiento para investigar el comportamiento de \(\widetilde{\beta}_0\). Para ello, considere que \(\widetilde{\beta}_0 = \bar{y} - \widetilde{\beta}_1\bar{x}\), por lo tanto:

\[E(\widetilde{\beta}_0)=E(\bar{y} - \widetilde{\beta}_1\bar{x})=E(\bar{y})-\bar{x}E(\widetilde{\beta}_1)=E(\sum_{i=1}^{n}\frac{y_i}{n})-\bar{x}\beta_1=\frac{1}{n}\sum_{i=1}^{n}E(y_i)-\bar{x}\beta_1\\=\frac{1}{n}\sum_{i=1}^{n}(\beta_0 + \beta_1x_i)-\bar{x}\beta_1=\frac{1}{n}\left(\sum_{i=1}^{n}\beta_0+\beta_1\sum_{i=1}^{n}x_i\right)-\bar{x}\beta_1\\=\frac{1}{n}(n\beta_0 + \beta_1n\bar{x})-\bar{x}\beta_1=\beta_0\]

De nuevo \(\widetilde{\beta}_0\) es un estimador insesgado de \(\beta_0\)!!!

Para la varianza tenemos:

\[Var(\widetilde{\beta}_0)=Var(\bar{y} - \widetilde{\beta}_1\bar{x})=Var(\bar{y})+{\bar{x}}^2Var(\widetilde{\beta}_1)=\frac{\sigma²}{n} + {\bar{x}}^2\frac{\sigma²}{S_{xx}}=\sigma²\left(\frac{1}{n}+\frac{{\bar{x}}^2}{S_{xx}}\right)\]

Una última anotación a realizar es que según el teorema de Gauss-Markov, para el modelo de regresión lineal con los supuestos \(E(e)=0\) y \(Var(e)=\sigma²\) y por supuesto con errores no correlacionados, los estimadores por mínimos cuadrados son insesgados (tal como acabamos de demostrar), pero más importante que eso, tienen la varianza más pequeña que se pueda esperar dentro del universo de estimadores insesgados que se construyen como combinaciones lineales de los \(y_i\). Es por esta razón que con frecuencia se dice que los estimadores por mínimos cuadrados son los estimadores lineales insesgados óptimos, donde óptimos implica que son de varianza mínima.

6.2 Simulando el comportamiento estadístico de los estimadores

Y exploramemos como se comportan los valores de \(\widetilde{\beta}_0\) y \(\widetilde{\beta}_1\) a lo largo de diversas muestras de tamaño \(n\) obtenidas con el modelo teórico. El pseudo código es:

  1. Generar una función para extraer muestras de tamaño \(n\) bajo un modelo teórico dado.
  2. Generar una función para ajustar una recta a la muestra anterior usando los estimadores \(\widetilde{\beta}_0\) y \(\widetilde{\beta}_1\) dados en las ecuaciones anteriores.
  3. Dibujar algunas rectas ajustadas para un conjunto de 5 muestras obtenidas del modelo teórico con \(\sigma=4\).
  4. Dibujar algunas rectas ajustadas para un conjunto de 5 muestras obtenidas del modelo teórico con \(\sigma=10\).
  5. Concluir según los resultados.
# Definimos un modelo teórico
n <- 50
vector_x <- c(1:n) * 30/50
sigma <- 4
beta_0 <- 2
beta_1 <- 0.5

# Creamos la función generadora de muestras
generador_muestra <-  function(n,vector_x,sigma, beta_0,beta_1){
  error <- rnorm(n,mean = 0, sd = sigma)
  y_observada <- beta_0 + beta_1 * vector_x + error
  return( list(x = vector_x, y = y_observada) )
}

# Corroborando que la función quedó bien 👍
generador_muestra(10,c(1:10),sigma, beta_0,beta_1)
## $x
##  [1]  1  2  3  4  5  6  7  8  9 10
## 
## $y
##  [1]  6.762204  2.392461  2.900142  6.312325  4.274832  3.802731  4.068061
##  [8]  7.206680  4.590557 -2.538485
# Creamos la función que estima una recta sobre una muestra dada
estimador_betas <- function(x,y){
    model <- lm(y~x)
    coeficientes <- model$coefficients
    return(list(beta_0 = coeficientes[1],
                beta_1 = coeficientes[2],
                sigma = sigma(model)))  
}

# Ajustamos semilla para reproducibilidad y una función que genera un
# gráfico con 5 rectas estimadas

set.seed(123)

generar_grafico <- function(x,n,sigma,beta_0,beta_1,position_legend){
df_rectas <- data.frame(x=numeric(),
                        y_pred=numeric(),
                        numero_muestra = numeric())
df_true <- data.frame(x=x,
                      y_true = (beta_0 + beta_1 * x))
for (i in 1:5) {
  muestra <- generador_muestra(n,x,sigma, beta_0,beta_1)
  betas_estimados <-  estimador_betas(muestra$x, muestra$y)
  df <- data.frame(x = x,
                   y_pred = (betas_estimados$beta_0 + 
                          betas_estimados$beta_1 * x),
                   numero_muestra = rep(i,n))
  df_rectas <- rbind(df_rectas,df)
}

p <- ggplot() +
  geom_line(data = df_rectas, aes(x=x, y=y_pred,
                                  color = factor(numero_muestra)))+
  scale_color_manual(values = c("1" = "blue", "2" = "black",
                                "3" = "green", "4" = "violet",
                                "5" = "orange",
                                "y real" = "red"),
                     name = "Rectas estimadas",
                     labels = c(paste("Muestra_",seq(1,5)),
                                "Recta real"))+
  labs(x="Valor x",
       y="y predicho",
       title = "Rectas estimadas sobre muestras distintas",
       subtitle = paste0("Sigma=",sigma, " n=",n))+
  geom_point(data = df_true, aes(x=x, y=y_true, color="y real"))+
  theme(legend.position=position_legend)
  return(p)
}

# Invocamos la función con valores de sigma distintos (4 y 10) y
# presentamos un gráfico comparativo.

grid.arrange(generar_grafico(x=vector_x,n=50,sigma=4,beta_0,beta_1,
                             position_legend='right'),
             generar_grafico(x=vector_x,n=50,sigma=10,beta_0,beta_1,
                             position_legend='right'),
             nrow=2)

Conclusiones:

  • Las rectas estimadas se aproximan a la recta real, pero dicha cercanía dependerá del sigma (variabilidad inherente del proceso generador de datos)
  • Como las rectas se estan distanciando de la recta teórica conviene investigar qué lo provoca, esto es, cuánto se distancian los betas estimados de los betas teóricos.
  • Conviene también estudiar el comportamiento estadístico de los betas estimados, no sólo la distancia a la que quedan de los betas teóricos.

Observación:

Hasta ahora no contamos con un estimado para sigma, pero es muy importante tenerlo porque del sigma del proceso generador de datos, está dependiendo la varianza de cada uno de los betas estimados, tal como se aprecia en el gráfico, en el que las rectas se distancian más de la teórica conforme el \(\sigma\) desconocido del proceso crece.

\(\widetilde{\sigma}²_e=\sum_{i=1}^{n}\frac{(y_i-\tilde{y_i})²}{n-p}\)

Acá \(p\) representa la cantidad de parámetros en mi modelo,es decir que el denominador debe ser los grados de libertad de la suma de cuadrados de los residuales para garantizar la insesgadez del estimador \(\widetilde{\sigma}²_e\), entendiendo por esto la propiedad que tiene de que su promedio sobre todas las posibles muestras extraídas del proceso generador de datos sea igual a \(\sigma²\), esto es \(E(\widetilde{\sigma}²_e)=\sigma²\)

Ahora estamos en condiciones de investigar el comportamiento estadístico de \(\widetilde{\beta}_0\) y \(\widetilde{\beta}_1\), para ello procederemos bajo el siguiente pseudo-código.

  • Fijar un tamaño de muestra y fijar un proceso generador de datos (puede usar el proceso del modelo anterior con \(\sigma = 4\))
  • Invocar 10000 veces el generador de la muestra
  • Sobre cada muestra invocar el estimador de los betas
  • Graficar el valor de \(\widetilde{\beta}_0\) en un histograma guardando el vector de betas estimados.
  • Repetir los pasos anteriores para tamaños de muestra cada vez más grandes.
  • Repetir los pasos anteriores para \(\widetilde{\beta}_1\)
  • Concluir sobre cómo el tamaño de muestra afecta la magnitud del error de estimación de los betas.
# Parámetros para fijar el proceso generador de datos
sigma <- 4
beta_0 <- 2
beta_1 <- 0.5

# vectores donde guardaremos los resultados de los betas estimados para
# cada tamaño de muestra
vector_beta_0_est <- vector()
vector_beta_1_est <- vector()
vector_sigma_est <- vector()
tamanos <- vector()

# Ciclo de alimentación de los vectores

if(!file.exists("df_estimaciones.csv")){
j<-1
for (n in c(5,10,15,20,50,100,500,1000)){
  for (i in 1:10000){
    muestra <-generador_muestra(n=n,vector_x=seq(1,n),sigma=sigma,
                                beta_0=beta_0,beta_1=beta_1)
    betas_estimados <- estimador_betas(muestra$x, muestra$y)
    vector_beta_0_est[j] <- betas_estimados$beta_0
    vector_beta_1_est[j] <- betas_estimados$beta_1
    vector_sigma_est[j] <- betas_estimados$sigma
    tamanos[j] <- n
    j <- j+1
  }
}

# Construyo el dataframe de estimaciones y en el caso de sigma estimado
# cálculo el estadístico de prueba que se usa para pruebas de hipótesis
# sobre sigma.

df_estimaciones <- data.frame(
  beta_0_est = vector_beta_0_est,
  beta_1_est = vector_beta_1_est,
  sigma_est = vector_sigma_est,
  chi_2 = (tamanos -2)*(vector_sigma_est^2)/sigma^2,
  tamano = tamanos
)

write.csv(df_estimaciones,file="df_estimaciones.csv",
          row.names = FALSE, sep=",")
}else{
  df_estimaciones <- read.csv("df_estimaciones.csv", sep=",")
}

df_estimaciones %>% 
  head(10) %>% 
  ftable()

beta_0_est

beta_1_est

sigma_est

chi_2

tamano

6.290113

-0.2369054

2.344651

1.0307600

5

-5.046949

3.1809693

3.355561

2.1112102

5

7.742303

-2.0710337

2.551786

1.2209269

5

9.809527

-1.9463694

5.283334

5.2338043

5

5.139514

-0.9018973

4.907130

4.5149854

5

2.811184

-0.2295792

2.211230

0.9167882

5

5.651254

-0.9752699

4.911290

4.5226448

5

-2.933659

2.7576534

7.369526

10.1831075

5

-1.740665

1.7925884

2.577523

1.2456793

5

7.470557

-0.4456503

6.246944

7.3170574

5

p1 <- ggplot(data=df_estimaciones, aes(x=beta_0_est))+
  geom_histogram()+
  labs(y = "Frecuencia", x="Beta 0 estimado")+
  facet_wrap(vars(tamano), scales="free")

plot(p1)

p2 <- ggplot(data=df_estimaciones, aes(x=beta_1_est))+
  geom_histogram()+
  labs(y = "Frecuencia", x="Beta 1 estimado")+
  facet_wrap(vars(tamano), scales="free")

plot(p2)

p3 <- ggplot(data=df_estimaciones, aes(x=chi_2))+
  geom_histogram()+
  labs(y = "Frecuencia", x="(n-2)^MSE/Sigma^2")+
  facet_wrap(vars(tamano), scales="free")

plot(p3)

Conclusiones:

  • Los histogramas para \(\widetilde{\beta}_0\) y \(\widetilde{\beta}_1\) tienen forma acampanada lo que sugiere una relación con la distribución normal (no estamos diciendo que se conduzcan como una normal)

  • Las varianzas de \(\widetilde{\beta}_0\) y \(\widetilde{\beta}_1\) dependen de los tamaños de muestra (y del sigma del proceso generador de datos aunque esta parte queda para inspección por parte del alumno)

  • Los histogramas para \(\widetilde{\sigma}²_e\) no se dibujaron, en su lugar se dibujó un histograma para una transformación (o escalamiento) dada por:

\[\chi² = \frac{(n-2)*\widetilde{\sigma}²_e}{\sigma^2}=\frac{SSE}{\sigma^2}= \frac{1}{\sigma^2}\sum_{i=1}^{n}(y_i - \widetilde{y})²\]

Es decir, en lugar de dibujar los valores originales de \(\widetilde{\sigma}²_e\), dibujamos la suma escalada de los cuadrados de los residuales, donde el factor de escalamiento es precisamente la varianza de los errores del proceso generador de datos.

¿Por qué hicimos esto? Porque esta suma de cuadrados, es una suma de cuadrados de variables aleatorias normales (le queda al estudiante demostrar por qué), es decir que ésta será una distribución chi cuadrado. Para determinar sus grados de libertad tomaremos en cuenta la cantidad de cuadrados independientes en la suma, entendiendo por independientes los cuadrados que pueden variar libremente sin violar las ecuaciones que dieron lugar a \(\widetilde{\beta_0}\) y \(\widetilde{\beta_1}\)

  • Será necesario derivar expresiones para la varianza de \(\widetilde\beta_1\) y \(\widetilde\beta_2\). Para ello tome en cuenta que las ecuaciones que los definen son combinaciones lineales de las \(y_i\). Además las \(y_i\) se asumen independientes (puesto que los \(e_i\) lo son por suposición previa). También recuerde que las \(y_i\) son variables aleatorias normales ya que dependen de los errores \(e_i\) los cuales son variables aleatorias normales (de nuevo por suposición previa).

Acá es donde notamos la importancia del supuesto de normalidad de los errores, y el supuesto de independencia, pues estos son la garantía para proponer los test de hipótesis que vamos a ver a continuación.

6.2 Prueba de hipótesis para \(\beta_1\)

Considere el siguiente estadístico de prueba:

\[t_0 = \frac{\widetilde\beta_0 - \beta_0}{\sqrt{Var(\widetilde\beta_0)}}\] Dada la normalidad de \(\widetilde{\beta_0}\) y el hecho de que \(\widetilde{\beta_0}\) es un estimador insesgado de \(\beta_0\) podemos decir que \(\t_0\) se distribuye como una variable normal estándar. Pero como obviamente tampoco conocemos el verdadero valor de la varianza de \(\beta_0\), en su lugar tendremos que usar un estimado del mismo \(\widetilde{Var}(\widetilde{\beta_0})\) dado por:

con \(n-2\) grados de libertad.

Entonces, podemos finalmente decir que tenemos instrumentos analíticos para proponer prueba de hipótesis para \(\beta_0\), \(\beta_1\), y \(\sigma\)

Para \(\beta_0\) podemos usar el estadístico de prueba dado por:

\(t_0 = \frac{\widetilde\beta_0 - \beta_{00}}{se(\widetilde\beta_0)}\)

con \(se(\widetilde\beta_0)=\sqrt{(MSE*\left(\frac{1}{n}+\frac{\bar{x}²}{S_{xx}}\right)}\)

Así que proponemos la prueba para \(\beta_0\):

\(H_0: \beta_0 = \beta_{00}\) vs \(H_a: \beta_0\ne\beta_{00}\) donde se rechaza \(H_0\) si \(|t_0|>t_{\frac{\alpha}{2}, n-2}\)

\(t_1 = \frac{\widetilde\beta_1 - \beta_1}{se(\widetilde\beta_1)}\)

Para \(\beta_1\) podemos usar el estadístico de prueba dado por:

\(t_0 = \frac{\widetilde\beta_1 - \beta_{10}}{se(\widetilde\beta_1)}\)

con \(se(\widetilde\beta_1)=\sqrt{\frac{MSE}{S_{xx}}}\)

Así que proponemos la prueba para \(\beta_1\):

\(H_0: \beta_1 = \beta_{10}\) vs \(H_a: \beta_0\ne\beta_{10}\) donde se rechaza \(H_0\) si \(|t_0|>t_{\frac{\alpha}{2}, n-2}\)

Para \(\sigma²\) podemos usar el estadístico de prueba dado por:

Variable aleatoria: Una variable aleatoria X es una función de \[ f : \Omega \rightarrow \mathbb{R}\] con \(\Omega\) el conjunto de resultados posibles de un experimento aleatorio. En otras palabras un mapeo desde el espacio de descripciones muestrales a los números reales con el fin de asociar a los resultados del experimento un valor numérico con el cual habilitar operaciones matemáticas.

[Proceso estocástico]{style = “color: red;”}: Sea \(\mathcal{T}\) un conjunto arbitrario. Un proceso estocástico es una familia \(\{X(t), \text{ t} \in \tau \}\), tal que, para cada \(t \in \mathcal{T}, X(t)\) es una variable aleatoria. Es decir, es una secuencia de variables aleatorias ordenadas por un índice, típicamente tomado del conjunto de los números enteros, o del conjunto de los números reales.

[Proceso estocástico en tiempo discreto]{style = “color: red;}: Es un proceso estocástico con un índice temporal discreto. Significando que \(\tau \in \mathbb{Z}\)

[Especificación de un proceso estocástico]{style = “color: red;”}: Sea \(t_1, t_2, t_3,...,t_n\) elementos cualquiera de \(\mathcal{T}\) y consideremos

\[F(x_1, x_2, x_3, ..., x_n; t_1, t_2, t_3, ..., t_n)=P\{X(t_1)<=x_1, ..., X(t_n)<=x_n\}\] [Definición de un proceso estocástico estacionario]{style = “color: red;”}: Un proceso estocástico \(\{X(t), \text{ t} \in \mathcal{T} \}\) se dice estrictamente estacionario si todas las distribuciones finito dimensionales tal como se definen en la ecuación anterior permanecen siendo las mismas sobre traslaciones en el tiempo, es decir:

\[F(x_1, x_2, x_3, ..., x_n; t_1+\tau, t_2+\tau, t_3+\tau, ..., t_n+\tau)=F(x_1, x_2, x_3, ..., x_n; t_1, t_2, t_3, ..., t_n)\] Algunas consecuencias de lo anterior son:

  • \(E\{X(t)\} = \mu(t) = \mu\)
  • \(Var\{X(t)\} = \sigma²(t) = \sigma²\)
  • \(\gamma(t_1, t_2) = \gamma(t_1-t_2, 0) = Cov\{X(t_1-t_2), X(0)\}\) o equivalentemente \(\gamma{\tau} = Cov\{X(t), X(t+\tau)\} = Cov\{X(0), X(\tau)\}\) para t, \(\tau \in \mathcal{T}\)

Reflexiones en lenguaje informal:

  • Un proceso estocástico quedará suficientemente definido si definimos la distribución de probabilidad n-dimensional conjunta para cada \(n>=1\), es decir, tomando \(n\) variables del proceso estocástico a la vez, sería necesario especificar la relación estadística de esas n variables a la vez para cada eleción de n. En nuestro contexto implica especificar cómo se espera que el precio de hoy se relacione con el precio de ayer, o cómo se espera que el precio de hoy se relacione con el de ayer y además con el de antier. Y así sucesivamente para cualquier grupo finito de precios que decida tomar a consideración.

  • Un proceso estocástico estacionario se puede interpretar como aquel cuyas propiedades estadísticas permanencen constantes a lo largo del tiempo y sólo cambiaran en función del grupo de \(n\) variables que decida obervar conjuntamente, pero no de la posición de esas n variables en el tiempo. Por ejemplo, suponiendo que n=1, esto implicaría que las propiedades estadísticas de cada variable tomada individualmente no se van alteradas conforme el tiempo avance, en particular cada variable será estable en media y en varianza. Lo que significa que puedo reunir las observaciones individuales en cada momento en el tiempo para estimar la media constante y la varianza de toda la secuencia, en la media en que en cada punto temporal estoy asumiendo que la media y la varianza es igual. Esto siempre se cumplirá? Por supuesto que no, es una suposición que imponemos a los datos para poder hacer inferencia con un único set de valores históricos observados.

  • Si la condición de estacionariedad estricta falla no hay nada que se pueda hacer en el estado actual del avance de la ciencia estadística.🤔 La única alternativa es relajar un poco las condiciones para exigencia de estacionariedad, cómo por ejemplo exigir estacionariedad débil consistente en pedir sólo homogeneidad en media, en varianza y en covarianza de segundo orden a la serie de tiempo bajo estudio. Exigir estacionariedad estricta es algo muy complicado de garantizar en la naturaleza.

  • Además es más fácil diseñar pruebas de estacionariedad débil.

Apliquemos una para la serie de tiempo de residuales

# Prueba de Dickey-FUller Aumentada
resultado_adf <- adf.test(mod1$residuals)
print(resultado_adf)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  mod1$residuals
## Dickey-Fuller = -3.6959, Lag order = 5, p-value = 0.02729
## alternative hypothesis: stationary
# Prueba de White
resultado_white <- bptest(mod1)
# Muestra el resultado de la prueba
print(resultado_white)
## 
##  studentized Breusch-Pagan test
## 
## data:  mod1
## BP = 1.741, df = 1, p-value = 0.187

En este caso debemos interpretar valores p<0.05 de la prueba de Dickey como un apoyo a la hipótesis alternativa de estacionariedad. Es decir, podemos considerar los residuales del modelo como estacionarios, porque un modelo autoregresivo con 5 rezagos no logró encontrar evidencia estadística de autocorrelaciones del error conducentes a no estacionariedad.

Por otro lado, podemos interpretar valores p<0.05 de la prueba de White como un apoyo a la hipótesis de homocedasticidad, es decir, podemos considerar los residuales del modelo como constantes en varianza a lo largo del tiempo.

Varianza de los residuales del modelo

Obteniendo los coeficientes.

mod1$coefficients
##  (Intercept) precio_dolar 
## 1436.6790256    0.9449775

Valores ajustados - Fitted values.

mod1$fitted.values %>% 
  head(20)
##        1        2        3        4        5        6        7        8 
## 3130.297 3123.597 3143.631 3139.479 3143.000 3157.035 3131.911 3109.503 
##        9       10       11       12       13       14       15       16 
## 3129.884 3148.096 3165.846 3183.871 3242.079 3230.649 3237.020 3248.871 
##       17       18       19       20 
## 3218.826 3253.241 3264.442 3291.015

Residuales

mod1$residuals %>% 
  head(20)
##          1          2          3          4          5          6          7 
##  340.70260  335.90285  399.86948  652.52072  518.25021  423.46463 -185.91095 
##          8          9         10         11         12         13         14 
## -333.25250 -398.13417  164.70412 -127.34649 -238.37054 -245.07936  287.60071 
##         15         16         17         18         19         20 
##  581.57961  273.37874  246.42423  -73.04054 -439.44215 -171.26527
# Los valores ajustados y los residuales también se pueden recuperar usando las funciones fitted( ) y residuals( ). Consulte la ayuda de estas funciones para conocer otros detalles.

Calcular la varianza de los residuales.

# En esta sección el estudiante desarrollará los cálculos.

Porcentaje de variabilidad no explicada por el modelo.

# En esta sección el estudiante desarrollará los cálculos.

Diagrama de dispersión con los puntos originales

Creación de la columna de predicción

En este chunk trabajamos con el modo de escritura Tidyverse, aunque se muestra cómo sería con R base.

#hass_dolar$predicciones <- predict(mod1)

hass_dolar <- hass_dolar %>% 
  mutate(predicciones = predict(mod1))
hass_dolar %>% 
  ggplot(aes(x = precio_dolar, y = precio_aguacate_kg)) +
  geom_smooth(method = "lm", se = FALSE, color="lightblue") +
  geom_segment(aes(xend=precio_dolar, yend=predicciones),
               col='red', lty='dashed') +
  geom_point() +
  geom_point(aes(y=predicciones), col='red') +
  theme_light()

#Con R Base
# ggplot(datos, aes(x=Edad, y=Resistencia)) +
#   geom_smooth(method="lm", se=FALSE, color="lightgrey") +
#   geom_segment(aes(xend=Edad, yend=predicciones), col='red', lty='dashed') +
#   geom_point() +
#   geom_point(aes(y=predicciones), col='red') +
#   theme_light()

REVISIÓN BÁSICA DE CONCEPTOS - Simulación.

Vamos a simular un modelo de regresión cuya especificación funcional es la siguiente:

\[ y = \beta_0 + \beta_1 * x + e \]

Con \(e \text{~} N(0,\sigma^2)\)

Como se puede notar, todo modelo teorico se compone de dos términos:

  1. El valor E(y|x)= E_y_x (valor esperado de y dado x)
  2. El componente aleatorio también llamado error.

Aunque en el modelo anterior hemos elegido una estructura lineal para representar E(y|x) en realidad podemos elegir cualquier otra estructura alternativa, siempre que respetemos la linealidad en los betas.

Dibujamos un gráfico de dispersión

DISCUTIENDO SOBRE LA LINEALIDAD DE LOS BETAS

En esta sección vamos a aclarar el asunto de que los BETAS sean lineales.

La linealidad en este contexto hace referencia a que no tengan potencias y que los coeficientes sean independendientes entre ellos.

Supongamos un modelo con coeficientes repetidos.

\[E(y|x) = \beta_0 + \alpha * x + \alpha * x^2 \]

Debemos factorizarlos para evitar estimarlos por separado y producir inconsistencias,esto debido a que ambos alpha (coeficientes) son dependientes.Es decir, reescribimos:

\[E(y|x) = \beta_0 + \alpha (x + x^2) \]

Donde el componente x + x² será el vector de valores de la variable predictora, que puede pensarse como una variable nueva:

\[ \tilde{x} = x + x^2\]

Y por lo tanto pensar al modelo como:

\[ E(y|x) = \beta_0 + \alpha * \tilde{x}\]

Esta transformación garantiza la creación de un nuevo modelo lineal respecto a \(\beta_0\) y \(\alpha\).

Es importante anotar que el dataframe que alimentará el modelo deberá tener computada la variable auxiliar \(\tilde{x}\), la cual debe ser pasada a la función encargada de estimar parámetros.

Esto corresponde a crear una nueva variable con los cálculos mencionados \(\tilde{x} = x + x^2\). Sin embargo, para graficar es conveniente usar la variable original.

ESTIMACIÓN DE PARÁMETROS

Por máxima verosimilitud.

Por mínimos cuadrados.

Modelo lineal

Unión en un dataframe

Gráficando los valores

Diferencia entre esperado y el predicho (Error)

Construir función

Usando la función

Creando vector BETAS

Histograma comparando distribuciones.